diff --git a/.gitignore b/.gitignore
index 3226ba4..182aa47 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,3 +7,6 @@ result
 # Python
 __pycache__
 *.py[cod]
+
+# Model build output
+/build
diff --git a/README.md b/README.md
index 282f834..ae83ed0 100644
--- a/README.md
+++ b/README.md
@@ -15,3 +15,10 @@ and this should succeed
 python3 -c "import nhf"
 ```
 
+## Testing
+
+Run all tests with
+``` sh
+python3 -m unittest
+```
+
diff --git a/nhf/__init__.py b/nhf/__init__.py
index e69de29..b582e6a 100644
--- a/nhf/__init__.py
+++ b/nhf/__init__.py
@@ -0,0 +1 @@
+from nhf.materials import Material, Role
diff --git a/nhf/build.py b/nhf/build.py
new file mode 100644
index 0000000..4435508
--- /dev/null
+++ b/nhf/build.py
@@ -0,0 +1,148 @@
+"""
+The NHF build system
+
+Usage: For any parametric assembly, inherit the `Model` class, and mark the
+output objects with the `@target` decorator. Each marked function should only
+take `self` as an argument.
+```python
+class BuildScaffold(Model):
+
+    @target(name="obj1")
+    def o1(self):
+        return Cq.Solid.makeBox(10, 10, 10)
+
+    def o2(self):
+        return Cq.Solid.makeCylinder(10, 20)
+```
+"""
+from enum import Enum
+from pathlib import Path
+from typing import Union
+from functools import wraps
+from colorama import Fore, Style
+import cadquery as Cq
+
+class TargetKind(Enum):
+
+    STL = "stl",
+    DXF = "dxf",
+
+    def __init__(self, ext: str):
+        self.ext = ext
+
+class Target:
+
+    def __init__(self,
+                 method,
+                 name: str,
+                 kind: TargetKind = TargetKind.STL,
+                 **kwargs):
+        self._method = method
+        self.name = name
+        self.kind = kind
+        self.kwargs = kwargs
+    def __str__(self):
+        return f"<target {self.name}.{self.kind.ext} {self._method}>"
+    def __call__(self, obj, *args, **kwargs):
+        """
+        Raw call function which passes arguments directly to `_method`
+        """
+        return self._method(obj, *args, **kwargs)
+
+    @property
+    def file_name(self):
+        """
+        Output file name
+        """
+        return f"{self.name}.{self.kind.ext}"
+
+    def write_to(self, obj, path: str):
+        x = self._method(obj)
+        if self.kind == TargetKind.STL:
+            assert isinstance(x, Union[
+                Cq.Workplane, Cq.Shape, Cq.Compound, Cq.Assembly])
+            if isinstance(x, Cq.Workplane):
+                x = x.val()
+            if isinstance(x, Cq.Assembly):
+                x = x.toCompound()
+            x.exportStl(path, **self.kwargs)
+        elif self.kind == TargetKind.DXF:
+            assert isinstance(x, Cq.Workplane)
+            Cq.exporters.exportDXF(x, path, **self.kwargs)
+        else:
+            assert False, f"Invalid kind: {self.kind}"
+
+    @classmethod
+    def methods(cls, subject):
+        """
+        List of all methods of a class or objects annotated with this decorator.
+        """
+        def g():
+            for name in dir(subject):
+                if name == 'target_names':
+                    continue
+                method = getattr(subject, name)
+                if hasattr(method, '_target'):
+                    yield method._target
+        return {method.name: method for method in g()}
+
+
+def target(name, **deco_kwargs):
+    """
+    Decorator for annotating a build output
+    """
+    def f(method):
+        @wraps(method)
+        def wrapper(self, *args, **kwargs):
+            return method(self, *args, **kwargs)
+        wrapper._target = Target(method, name, **deco_kwargs)
+        return wrapper
+    return f
+
+
+class Model:
+    """
+    Base class for a parametric assembly
+    """
+    def __init__(self, name: str):
+        self.name = name
+
+    @property
+    def target_names(self) -> list[str]:
+        """
+        List of decorated target functions
+        """
+        return list(Target.methods(self).keys())
+
+    def check_all(self) -> int:
+        """
+        Builds all targets but do not output them
+        """
+        total = 0
+        for t in Target.methods(self).values():
+            t(self)
+            total += 1
+        return total
+
+    def build_all(self, output_dir: Union[Path, str] = "build", verbose=1):
+        """
+        Build all targets in this model and write the results to file
+        """
+        output_dir = Path(output_dir)
+        for t in Target.methods(self).values():
+            output_file = output_dir / self.name / t.file_name
+            if output_file.is_file():
+                if verbose >= 1:
+                    print(f"{Fore.GREEN}Skipping{Style.RESET_ALL} {output_file}")
+                continue
+            output_file.parent.mkdir(exist_ok=True, parents=True)
+
+            if verbose >= 1:
+                print(f"{Fore.BLUE}Building{Style.RESET_ALL} {output_file}")
+
+            try:
+                t.write_to(self, str(output_file))
+                if verbose >= 1:
+                    print(f"{Fore.GREEN}Built{Style.RESET_ALL} {output_file}")
+            except Exception as e:
+                print(f"{Fore.RED}Failed to build{Style.RESET_ALL} {output_file}: {e}")
diff --git a/nhf/checks.py b/nhf/checks.py
new file mode 100644
index 0000000..468cc8e
--- /dev/null
+++ b/nhf/checks.py
@@ -0,0 +1,6 @@
+import cadquery as Cq
+
+def binary_intersection(a: Cq.Assembly) -> Cq.Shape:
+    objs = [s.toCompound() for _, s in a.traverse() if isinstance(s, Cq.Assembly)]
+    obj1, obj2 = objs[:2]
+    return obj1.intersect(obj2)
diff --git a/nhf/diag.py b/nhf/diag.py
new file mode 100644
index 0000000..1648383
--- /dev/null
+++ b/nhf/diag.py
@@ -0,0 +1,42 @@
+import cadquery as Cq
+
+def tidy_repr(obj):
+    """Shortens a default repr string"""
+    return repr(obj).split(".")[-1].rstrip(">")
+
+
+def _ctx_str(self):
+    return (
+        tidy_repr(self)
+        + ":\n"
+        + f"    pendingWires: {self.pendingWires}\n"
+        + f"    pendingEdges: {self.pendingEdges}\n"
+        + f"    tags: {self.tags}"
+    )
+
+
+Cq.cq.CQContext.__str__ = _ctx_str
+
+
+def _plane_str(self):
+    return (
+        tidy_repr(self)
+        + ":\n"
+        + f"    origin: {self.origin.toTuple()}\n"
+        + f"    z direction: {self.zDir.toTuple()}"
+    )
+
+
+Cq.occ_impl.geom.Plane.__str__ = _plane_str
+
+
+def _wp_str(self):
+    out = tidy_repr(self) + ":\n"
+    out += f"  parent: {tidy_repr(self.parent)}\n" if self.parent else "  no parent\n"
+    out += f"  plane: {self.plane}\n"
+    out += f"  objects: {self.objects}\n"
+    out += f"  modelling context: {self.ctx}"
+    return out
+
+
+Cq.Workplane.__str__ = _wp_str
diff --git a/nhf/materials.py b/nhf/materials.py
new file mode 100644
index 0000000..c0adbde
--- /dev/null
+++ b/nhf/materials.py
@@ -0,0 +1,38 @@
+"""
+A catalog of material properties
+"""
+from enum import Enum
+import cadquery as Cq
+
+def _color(name: str, alpha: float) -> Cq.Color:
+    r, g, b, _ = Cq.Color(name).toTuple()
+    return Cq.Color(r, g, b, alpha)
+
+class Role(Enum):
+    """
+    Describes the role of a part
+    """
+
+    # Parent and child components in a load bearing joint
+    PARENT = _color('blue4', 0.6)
+    CHILD = _color('darkorange2', 0.6)
+    DAMPING = _color('springgreen', 0.5)
+    STRUCTURE = _color('gray', 0.4)
+    DECORATION = _color('lightseagreen', 0.4)
+    ELECTRONIC = _color('mediumorchid', 0.5)
+
+    def __init__(self, color: Cq.Color):
+        self.color = color
+
+class Material(Enum):
+    """
+    A catalog of common material properties
+    """
+
+    WOOD_BIRCH = 0.8, _color('bisque', 0.9)
+    PLASTIC_PLA = 0.5, _color('azure3', 0.6)
+    ACRYLIC_BLACK = 0.5, _color('gray50', 0.6)
+
+    def __init__(self, density: float, color: Cq.Color):
+        self.density = density
+        self.color = color
diff --git a/nhf/parts/__init__.py b/nhf/parts/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/nhf/parts/handle.py b/nhf/parts/handle.py
new file mode 100644
index 0000000..c6845ef
--- /dev/null
+++ b/nhf/parts/handle.py
@@ -0,0 +1,252 @@
+"""
+This schematics file contains all designs related to tool handles
+"""
+from dataclasses import dataclass
+import cadquery as Cq
+import nhf.parts.metric_threads as metric_threads
+
+@dataclass
+class Handle:
+    """
+    Characteristic of a tool handle
+
+    This assumes the handle segment material does not have threads. Each segment
+    attaches to two insertions, which have threads on the inside. A connector
+    has threads on the outside and joints two insertions.
+
+    Note that all the radial sizes are diameters (in mm).
+    """
+
+    # Outer and inner radius for the handle usually come in standard sizes
+    diam: float = 38
+    diam_inner: float = 33
+
+    # Major diameter of the internal threads, following ISO metric screw thread
+    # standard. This determines the wall thickness of the insertion.
+    diam_threading: float = 27.0
+
+    thread_pitch: float = 3.0
+
+    # Internal cavity diameter. This determines the wall thickness of the connector
+    diam_connector_internal: float = 18.0
+
+    # If set to true, do not generate threads
+    simplify_geometry: bool = True
+
+    # Length for the rim on the female connector
+    rim_length: float = 5
+
+    insertion_length: float = 30
+
+    # Amount by which the connector goes into the segment
+    connector_length: float = 60
+
+    def __post_init__(self):
+        assert self.diam > self.diam_inner, "Material thickness cannot be <= 0"
+        assert self.diam_inner > self.diam_insertion_internal, "Threading radius is too big"
+        assert self.diam_insertion_internal > self.diam_connector_external
+        assert self.diam_connector_external > self.diam_connector_internal, "Internal diameter is too large"
+        assert self.insertion_length > self.rim_length
+
+    @property
+    def diam_insertion_internal(self):
+        r = metric_threads.metric_thread_major_radius(
+            self.diam_threading,
+            self.thread_pitch,
+            internal=True)
+        return r * 2
+
+    @property
+    def diam_connector_external(self):
+        r = metric_threads.metric_thread_minor_radius(
+            self.diam_threading,
+            self.thread_pitch)
+        return r * 2
+
+    def segment(self, length: float):
+        result = (
+            Cq.Workplane()
+            .cylinder(
+                radius=self.diam / 2,
+                height=length)
+        )
+        result.faces("<Z").tag("mate1")
+        result.faces(">Z").tag("mate2")
+        return result
+
+    def _external_thread(self, length=None):
+        if length is None:
+            length = self.insertion_length
+        return metric_threads.external_metric_thread(
+            self.diam_threading,
+            self.thread_pitch,
+            length,
+            top_lead_in=True)
+    def _internal_thread(self):
+        return metric_threads.internal_metric_thread(
+            self.diam_threading,
+            self.thread_pitch,
+            self.insertion_length)
+
+    def insertion(self, holes=[]):
+        """
+        This type of joint is used to connect two handlebar pieces. Each handlebar
+        piece is a tube which cannot be machined, so the joint connects to the
+        handle by glue.
+
+        Tags:
+        * lip: Co-planar Mates to the rod
+        * mate: Mates to the connector
+
+        WARNING: A tolerance lower than the defualt (maybe 5e-4) is required for
+        STL export.
+
+        Set `holes` to the heights for drilling holes into the model for resin
+        to flow out.
+        """
+        result = (
+            Cq.Workplane('XY')
+            .cylinder(
+                radius=self.diam_inner / 2,
+                height=self.insertion_length - self.rim_length,
+                centered=[True, True, False])
+        )
+        result.faces(">Z").tag("rim")
+        if self.rim_length > 0:
+            result = (
+                result.faces(">Z")
+                .workplane()
+                .circle(self.diam / 2)
+                .extrude(self.rim_length)
+                .faces(">Z")
+                .hole(self.diam_insertion_internal)
+            )
+        result.faces(">Z").tag("mate")
+        if not self.simplify_geometry:
+            thread = self._internal_thread().val()
+            result = result.union(thread)
+        for h in holes:
+            cyl = Cq.Solid.makeCylinder(
+                radius=2,
+                height=self.diam * 2,
+                pnt=(-self.diam, 0, h),
+                dir=(1, 0, 0))
+            result = result.cut(cyl)
+        return result
+
+    def connector(self, solid: bool = True):
+        """
+        Tags:
+        * mate{1,2}: Mates to the connector
+
+        WARNING: A tolerance lower than the defualt (maybe 2e-4) is required for
+        STL export.
+        """
+        result = (
+            Cq.Workplane('XY')
+            .cylinder(
+                radius=self.diam / 2,
+                height=self.connector_length,
+            )
+        )
+        for (tag, selector) in [("mate1", "<Z"), ("mate2", ">Z")]:
+            result.faces(selector).tag(tag)
+            result = (
+                result
+                .faces(selector)
+                .workplane()
+                .circle(self.diam_connector_external / 2)
+                .extrude(self.insertion_length)
+            )
+        if not solid:
+            result = result.faces(">Z").hole(self.diam_connector_internal)
+        if not self.simplify_geometry:
+            thread = self._external_thread().val()
+            result = (
+                result
+                .union(
+                    thread
+                    .located(Cq.Location((0, 0, self.connector_length / 2))))
+                .union(
+                    thread
+                    .rotate((0,0,0), (1,0,0), angleDegrees=180)
+                    .located(Cq.Location((0, 0, -self.connector_length / 2))))
+            )
+        return result
+
+    def one_side_connector(self, height=None):
+        if height is None:
+            height = self.rim_length
+        result = (
+            Cq.Workplane('XY')
+            .cylinder(
+                radius=self.diam / 2,
+                height=height,
+                centered=(True, True, False)
+            )
+        )
+        result.faces(">Z").tag("mate")
+        result.faces("<Z").tag("base")
+        result = (
+            result
+            .faces(">Z")
+            .workplane()
+            .circle(self.diam_connector_external / 2)
+            .extrude(self.insertion_length)
+        )
+        if not self.simplify_geometry:
+            thread = self._external_thread().val()
+            result = (
+                result
+                .union(
+                    thread
+                    .located(Cq.Location((0, 0, height))))
+            )
+        return result
+
+    def threaded_core(self, length):
+        """
+        Generates a threaded core for unioning with other components
+        """
+        result = (
+            Cq.Workplane('XY')
+            .cylinder(
+                radius=self.diam_connector_external / 2,
+                height=length,
+                centered=(True, True, False),
+            )
+        )
+        result.faces(">Z").tag("mate")
+        result.faces("<Z").tag("base")
+        if not self.simplify_geometry:
+            thread = self._external_thread(length=length).val()
+            result = (
+                result
+                .union(thread)
+            )
+        return result
+
+    def connector_insertion_assembly(self):
+        connector_color = Cq.Color(0.8,0.8,0.5,0.3)
+        insertion_color = Cq.Color(0.7,0.7,0.7,0.3)
+        result = (
+            Cq.Assembly()
+            .add(self.connector(), name="c", color=connector_color)
+            .add(self.insertion(), name="i1", color=insertion_color)
+            .add(self.insertion(), name="i2", color=insertion_color)
+            .constrain("c?mate1", "i1?mate", "Plane")
+            .constrain("c?mate2", "i2?mate", "Plane")
+            .solve()
+        )
+        return result
+    def connector_one_side_insertion_assembly(self):
+        connector_color = Cq.Color(0.8,0.8,0.5,0.3)
+        insertion_color = Cq.Color(0.7,0.7,0.7,0.3)
+        result = (
+            Cq.Assembly()
+            .add(self.insertion(), name="i", color=connector_color)
+            .add(self.one_side_connector(), name="c", color=insertion_color)
+            .constrain("i?mate", "c?mate", "Plane")
+            .solve()
+        )
+        return result
diff --git a/nhf/parts/joints.py b/nhf/parts/joints.py
new file mode 100644
index 0000000..036a4f4
--- /dev/null
+++ b/nhf/parts/joints.py
@@ -0,0 +1,421 @@
+from dataclasses import dataclass
+import math
+import cadquery as Cq
+import nhf.parts.springs as springs
+from nhf import Role
+
+@dataclass
+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
+
+    def __post_init__(self):
+        # Ensures tangent doesn't blow up
+        assert self.n_tooth >= 5
+        assert self.radius > self.radius_inner
+
+    @property
+    def tooth_angle(self):
+        return 360 / self.n_tooth
+
+    @property
+    def total_height(self):
+        return self.base_height + self.tooth_height
+
+
+    def generate(self, is_mated=False, tol=0.01):
+        """
+        is_mated: If set to true, rotate the teeth so they line up at 0 degrees.
+
+        FIXME: Mate is not exact when number of tooth is low
+        """
+        phi = math.radians(self.tooth_angle)
+        alpha = 2 * math.atan(self.radius / self.tooth_height * math.tan(phi/2))
+        #alpha = math.atan(self.radius * math.radians(180 / self.n_tooth) / self.tooth_height)
+        gamma = math.radians(90 / self.n_tooth)
+        # Tooth half height
+        l = self.radius * math.cos(gamma)
+        a = self.radius * math.sin(gamma)
+        t = a / math.tan(alpha / 2)
+        beta = math.asin(t / l)
+        dx = self.tooth_height * math.tan(alpha / 2)
+        profile = (
+            Cq.Workplane('YZ')
+            .polyline([
+                (0, 0),
+                (dx, self.tooth_height),
+                (-dx, self.tooth_height),
+            ])
+            .close()
+            .extrude(-self.radius)
+            .val()
+            .rotate((0, 0, 0), (0, 1, 0), math.degrees(beta))
+            .moved(Cq.Location((0, 0, self.base_height)))
+        )
+        core = Cq.Solid.makeCylinder(
+            radius=self.radius_inner,
+            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("bore")
+            .cut(core)
+            .polarArray(
+                radius=self.radius,
+                startAngle=angle_offset,
+                angle=360,
+                count=self.n_tooth)
+            .cutEach(
+                lambda loc: profile.moved(loc),
+            )
+        )
+        (
+            result
+            .polyline([
+                (0, 0, self.base_height),
+                (0, 0, self.base_height + self.tooth_height)
+            ], forConstruction=True)
+            .tag("mate")
+        )
+        (
+            result
+            .polyline([(0, 0, 0), (1, 0, 0)], forConstruction=True)
+            .tag("directrix")
+        )
+        return result
+
+    def assembly(self, offset: int = 1):
+        """
+        Generate an example assembly
+        """
+        tab = (
+            Cq.Workplane('XY')
+            .box(100, 10, 2, centered=False)
+        )
+        obj1 = (
+            self.generate()
+            .faces(tag="bore")
+            .cboreHole(
+                diameter=10,
+                cboreDiameter=20,
+                cboreDepth=3)
+            .union(tab)
+        )
+        obj2 = (
+            self.generate(is_mated=True)
+            .union(tab)
+        )
+        angle = offset * self.tooth_angle
+        result = (
+            Cq.Assembly()
+            .add(obj1, name="obj1", color=Role.PARENT.color)
+            .add(obj2, name="obj2", color=Role.CHILD.color)
+            .constrain("obj1", "Fixed")
+            .constrain("obj1?mate", "obj2?mate", "Plane")
+            .constrain("obj1?directrix", "obj2?directrix", "Axis", param=angle)
+            .solve()
+        )
+        return result
+
+def comma_joint(radius=30,
+                shaft_radius=10,
+                height=10,
+                flange=10,
+                flange_thickness=25,
+                n_serration=16,
+                serration_angle_offset=0,
+                serration_height=5,
+                serration_inner_radius=20,
+                serration_theta=2 * math.pi / 48,
+                serration_tilt=-30,
+                right_handed=False):
+    """
+    Produces a "o_" shaped joint, with serrations to accomodate a torsion spring
+    """
+    assert flange_thickness <= radius
+    flange_poly = [
+        (0, radius - flange_thickness),
+        (0, radius),
+        (flange + radius, radius),
+        (flange + radius, radius - flange_thickness)
+    ]
+    if right_handed:
+        flange_poly = [(x, -y) for x,y in flange_poly]
+    sketch = (
+        Cq.Sketch()
+        .circle(radius)
+        .polygon(flange_poly, mode='a')
+        .circle(shaft_radius, mode='s')
+    )
+    serration_poly = [
+        (0, 0), (radius, 0),
+        (radius, radius * math.tan(serration_theta))
+    ]
+    serration = (
+        Cq.Workplane('XY')
+        .sketch()
+        .polygon(serration_poly)
+        .circle(radius, mode='i')
+        .circle(serration_inner_radius, mode='s')
+        .finalize()
+        .extrude(serration_height)
+        .translate(Cq.Vector((-serration_inner_radius, 0, height)))
+        .rotate(
+            axisStartPoint=(0, 0, 0),
+            axisEndPoint=(0, 0, height),
+            angleDegrees=serration_tilt)
+        .val()
+    )
+    serrations = (
+        Cq.Workplane('XY')
+        .polarArray(radius=serration_inner_radius,
+                    startAngle=0+serration_angle_offset,
+                    angle=360+serration_angle_offset,
+                    count=n_serration)
+        .eachpoint(lambda loc: serration.located(loc))
+    )
+    result = (
+        Cq.Workplane()
+        .add(sketch)
+        .extrude(height)
+        .union(serrations)
+        .clean()
+    )
+
+    result.polyline([
+        (0, 0, height - serration_height),
+        (0, 0, height + serration_height)],
+        forConstruction=True).tag("serrated")
+    result.polyline([
+        (0, radius, 0),
+        (flange + radius, radius, 0)],
+        forConstruction=True).tag("tail")
+    result.faces('>X').tag("tail_end")
+    return result
+
+def comma_assembly():
+    joint1 = comma_joint()
+    joint2 = comma_joint()
+    spring = springs.torsion_spring()
+    result = (
+        Cq.Assembly()
+        .add(joint1, name="joint1", color=Cq.Color(0.8,0.8,0.5,0.3))
+        .add(joint2, name="joint2", color=Cq.Color(0.8,0.8,0.5,0.3))
+        .add(spring, name="spring", color=Cq.Color(0.5,0.5,0.5,1))
+        .constrain("joint1?serrated", "spring?bot", "Plane")
+        .constrain("joint2?serrated", "spring?top", "Plane")
+        .constrain("joint1?tail", "FixedAxis", (1, 0, 0))
+        .constrain("joint2?tail", "FixedAxis", (-1, 0, 0))
+        .solve()
+    )
+    return result
+
+@dataclass
+class TorsionJoint:
+    """
+    This jonit consists of a rider puck on a track puck. IT is best suited if
+    the radius has to be small and vertical space is abundant.
+    """
+
+    # Radius limit for rotating components
+    radius: float = 40
+    disk_height: float = 10
+
+    radius_spring: float = 15
+    radius_axle: float = 6
+
+    # Offset of the spring hole w.r.t. surface
+    spring_hole_depth: float = 4
+
+    # Also used for the height of the hole for the spring
+    spring_thickness: float = 2
+    spring_height: float = 15
+
+    spring_tail_length: float = 40
+
+    groove_radius_outer: float = 35
+    groove_radius_inner: float = 20
+    groove_depth: float = 5
+    rider_gap: float = 2
+    n_slots: float = 8
+
+    right_handed: bool = False
+
+    def __post_init__(self):
+        assert self.disk_height > self.spring_hole_depth
+        assert self.radius > self.groove_radius_outer
+        assert self.groove_radius_outer > self.groove_radius_inner
+        assert self.groove_radius_inner > self.radius_spring
+        assert self.spring_height > self.groove_depth, "Groove is too deep"
+        assert self.radius_spring > self.radius_axle
+
+    @property
+    def total_height(self):
+        return 2 * self.disk_height + self.spring_height
+
+    @property
+    def _radius_spring_internal(self):
+        return self.radius_spring - self.spring_thickness
+
+    def _slot_polygon(self, flip: bool=False):
+        r1 = self.radius_spring - self.spring_thickness
+        r2 = self.radius_spring
+        flip = flip != self.right_handed
+        if flip:
+            r1 = -r1
+            r2 = -r2
+        return [
+            (0, r2),
+            (self.spring_tail_length, r2),
+            (self.spring_tail_length, r1),
+            (0, r1),
+        ]
+    def _directrix(self, height, theta=0):
+        c, s = math.cos(theta), math.sin(theta)
+        r2 = self.radius_spring
+        l = self.spring_tail_length
+        if self.right_handed:
+            r2 = -r2
+        # This is (0, r2) and (l, r2) transformed by rotation matrix
+        # [[c, s], [-s, c]]
+        return [
+            (s * r2,         -s * l + c * r2, height),
+            (c * l + s * r2, -s * l + c * r2, height),
+        ]
+
+
+    def spring(self):
+        return springs.torsion_spring(
+            radius=self.radius_spring,
+            height=self.spring_height,
+            thickness=self.spring_thickness,
+            tail_length=self.spring_tail_length,
+        )
+
+    def track(self):
+        groove_profile = (
+            Cq.Sketch()
+            .circle(self.radius)
+            .circle(self.groove_radius_outer, mode='s')
+            .circle(self.groove_radius_inner, mode='a')
+            .circle(self.radius_spring, mode='s')
+        )
+        spring_hole_profile = (
+            Cq.Sketch()
+            .circle(self.radius)
+            .polygon(self._slot_polygon(flip=False), mode='s')
+            .circle(self.radius_spring, mode='s')
+        )
+        result = (
+            Cq.Workplane('XY')
+            .cylinder(
+                radius=self.radius,
+                height=self.disk_height,
+                centered=(True, True, False))
+            .faces('>Z')
+            .tag("spring")
+            .placeSketch(spring_hole_profile)
+            .extrude(self.spring_thickness)
+            # If the spring hole profile is not simply connected, this workplane
+            # will have to be created from the `spring-mate` face.
+            .faces('>Z')
+            .placeSketch(groove_profile)
+            .extrude(self.groove_depth)
+            .faces('>Z')
+            .hole(self.radius_axle * 2)
+        )
+        # Insert directrix`
+        result.polyline(self._directrix(self.disk_height),
+                        forConstruction=True).tag("directrix")
+        return result
+
+    def rider(self):
+        def slot(loc):
+            wire = Cq.Wire.makePolygon(self._slot_polygon(flip=False))
+            face = Cq.Face.makeFromWires(wire)
+            return face.located(loc)
+        wall_profile = (
+            Cq.Sketch()
+            .circle(self.radius, mode='a')
+            .circle(self.radius_spring, mode='s')
+            .parray(
+                r=0,
+                a1=0,
+                da=360,
+                n=self.n_slots)
+            .each(slot, mode='s')
+            #.circle(self._radius_wall, mode='a')
+        )
+        contact_profile = (
+            Cq.Sketch()
+            .circle(self.groove_radius_outer, mode='a')
+            .circle(self.groove_radius_inner, mode='s')
+            #.circle(self._radius_wall, mode='a')
+            .parray(
+                r=0,
+                a1=0,
+                da=360,
+                n=self.n_slots)
+            .each(slot, mode='s')
+        )
+        middle_height = self.spring_height - self.groove_depth - self.rider_gap
+        result = (
+            Cq.Workplane('XY')
+            .cylinder(
+                radius=self.radius,
+                height=self.disk_height,
+                centered=(True, True, False))
+            .faces('>Z')
+            .tag("spring")
+            .placeSketch(wall_profile)
+            .extrude(middle_height)
+            # The top face might not be in one piece.
+            #.faces('>Z')
+            .workplane(offset=middle_height)
+            .placeSketch(contact_profile)
+            .extrude(self.groove_depth + self.rider_gap)
+            .faces(tag="spring")
+            .circle(self._radius_spring_internal)
+            .extrude(self.spring_height)
+            .faces('>Z')
+            .hole(self.radius_axle * 2)
+        )
+        for i in range(self.n_slots):
+            theta = 2 * math.pi * i / self.n_slots
+            result.polyline(self._directrix(self.disk_height, theta),
+                            forConstruction=True).tag(f"directrix{i}")
+        return result
+
+    def rider_track_assembly(self):
+        rider = self.rider()
+        track = self.track()
+        spring = self.spring()
+        result = (
+            Cq.Assembly()
+            .add(spring, name="spring", color=Role.DAMPING.color)
+            .add(track, name="track", color=Role.PARENT.color)
+            .constrain("track?spring", "spring?top", "Plane")
+            .add(rider, name="rider", color=Role.CHILD.color)
+            .constrain("rider?spring", "spring?bot", "Plane")
+            .constrain("track?directrix", "spring?directrix_bot", "Axis")
+            .constrain("rider?directrix0", "spring?directrix_top", "Axis")
+            .solve()
+        )
+        return result
diff --git a/nhf/parts/metric_threads.py b/nhf/parts/metric_threads.py
new file mode 100644
index 0000000..3712fe0
--- /dev/null
+++ b/nhf/parts/metric_threads.py
@@ -0,0 +1,422 @@
+# Copyright (c) 2020-2024, Nerius Anthony Landys.  All rights reserved.
+#   neri-engineering 'at' protonmail.com
+#   https://svn.code.sf.net/p/nl10/code/cq-code/common/metric_threads.py
+# This file is public domain.  Use it for any purpose, including commercial
+# applications.  Attribution would be nice, but is not required.  There is no
+# warranty of any kind, including its correctness, usefulness, or safety.
+#
+# Simple code example to create meshing M3x0.5 threads:
+###############################################################################
+#
+#     male   = external_metric_thread(3.0, 0.5, 4.0, z_start= -0.85,
+#                                     top_lead_in=True)
+#
+#     # Please note that the female thread is meant for a hole which has
+#     # radius equal to metric_thread_major_radius(3.0, 0.5, internal=True),
+#     # which is in fact very slightly larger than a 3.0 diameter hole.
+#
+#     female = internal_metric_thread(3.0, 0.5, 1.5,
+#                                     bottom_chamfer=True, base_tube_od= 4.5)
+#
+###############################################################################
+# Left hand threads can be created by employing one of the "mirror" operations.
+# Thanks for taking the time to understand and use this code!
+
+import math
+import cadquery as cq
+
+###############################################################################
+# The functions which have names preceded by '__' are not meant to be called
+# externally; the remaining functions are written with the intention that they
+# will be called by external code.  The first section of code consists of
+# lightweight helper functions; the meat and potatoes of this library is last.
+###############################################################################
+
+# Return value is in degrees, and currently it's fixed at 30.  Essentially this
+# results in a typical 60 degree equilateral triangle cutting bit for threads.
+def metric_thread_angle():
+    return 30
+
+# Helper func. to make code more intuitive and succinct.  Degrees --> radians.
+def __deg2rad(degrees):
+    return degrees * math.pi / 180
+
+# In the absence of flat thread valley and flattened thread tip, returns the
+# amount by which the thread "triangle" protrudes outwards (radially) from base
+# cylinder in the case of external thread, or the amount by which the thread
+# "triangle" protrudes inwards from base tube in the case of internal thread.
+def metric_thread_perfect_height(pitch):
+    return pitch / (2 * math.tan(__deg2rad(metric_thread_angle())))
+
+# Up the radii of internal (female) thread in order to provide a little bit of
+# wiggle room around male thread.  Right now input parameter 'diameter' is
+# ignored.  This function is only used for internal/female threads.  Currently
+# there is no practical way to adjust the male/female thread clearance besides
+# to manually edit this function.  This design route was chosen for the sake of
+# code simplicity.
+def __metric_thread_internal_radius_increase(diameter, pitch):
+    return 0.1 * metric_thread_perfect_height(pitch)
+
+# Returns the major radius of thread, which is always the greater of the two.
+def metric_thread_major_radius(diameter, pitch, internal=False):
+    return (__metric_thread_internal_radius_increase(diameter, pitch) if
+            internal else 0.0) + (diameter / 2)
+
+# What portion of the total pitch is taken up by the angled thread section (and
+# not the squared off valley and tip).  The remaining portion (1 minus ratio)
+# will be divided equally between the flattened valley and flattened tip.
+def __metric_thread_effective_ratio():
+    return 0.7
+
+# Returns the minor radius of thread, which is always the lesser of the two.
+def metric_thread_minor_radius(diameter, pitch, internal=False):
+    return (metric_thread_major_radius(diameter, pitch, internal)
+            - (__metric_thread_effective_ratio() *
+               metric_thread_perfect_height(pitch)))
+
+# What the major radius would be if the cuts were perfectly triangular, without
+# flat spots in the valleys and without flattened tips.
+def metric_thread_perfect_major_radius(diameter, pitch, internal=False):
+    return (metric_thread_major_radius(diameter, pitch, internal)
+            + ((1.0 - __metric_thread_effective_ratio()) *
+               metric_thread_perfect_height(pitch) / 2))
+
+# What the minor radius would be if the cuts were perfectly triangular, without
+# flat spots in the valleys and without flattened tips.
+def metric_thread_perfect_minor_radius(diameter, pitch, internal=False):
+    return (metric_thread_perfect_major_radius(diameter, pitch, internal)
+            - metric_thread_perfect_height(pitch))
+
+# Returns the lead-in and/or chamfer distance along the z axis of rotation.
+# The lead-in/chamfer only depends on the pitch and is made with the same angle
+# as the thread, that being 30 degrees offset from radial.
+def metric_thread_lead_in(pitch, internal=False):
+    return (math.tan(__deg2rad(metric_thread_angle()))
+            * (metric_thread_major_radius(256.0, pitch, internal)
+               - metric_thread_minor_radius(256.0, pitch, internal)))
+
+# Returns the width of the flat spot in thread valley of a standard thread.
+# This is also equal to the width of the flat spot on thread tip, on a standard
+# thread.
+def metric_thread_relief(pitch):
+    return (1.0 - __metric_thread_effective_ratio()) * pitch / 2
+
+
+###############################################################################
+# A few words on modules external_metric_thread() and internal_metric_thread().
+# The parameter 'z_start' is added as a convenience in order to make the male
+# and female threads align perfectly.  When male and female threads are created
+# having the same diameter, pitch, and n_starts (usually 1), then so long as
+# they are not translated or rotated (or so long as they are subjected to the
+# same exact translation and rotation), they will intermesh perfectly,
+# regardless of the value of 'z_start' used on each.  This is in order that
+# assemblies be able to depict perfectly aligning threads.
+
+# Generates threads with base cylinder unless 'base_cylinder' is overridden.
+# Please note that 'use_epsilon' is activated by default, which causes a slight
+# budge in the minor radius, inwards, so that overlaps would be created with
+# inner cylinders.  (Does not affect thread profile outside of cylinder.)
+###############################################################################
+def external_metric_thread(diameter,  # Required parameter, e.g. 3.0 for M3x0.5
+                           pitch,     # Required parameter, e.g. 0.5 for M3x0.5
+                           length,    # Required parameter, e.g. 2.0
+                           z_start=0.0,
+                           n_starts=1,
+                           bottom_lead_in=False, # Lead-in is at same angle as
+                           top_lead_in   =False, # thread, namely 30 degrees.
+                           bottom_relief=False, # Add relief groove to start or
+                           top_relief   =False, # end of threads (shorten).
+                           force_outer_radius=-1.0, # Set close to diameter/2.
+                           use_epsilon=True, # For inner cylinder overlap.
+                           base_cylinder=True, # Whether to include base cyl.
+                           cyl_extend_bottom=-1.0,
+                           cyl_extend_top=-1.0,
+                           envelope=False): # Draw only envelope, don't cut.
+
+    cyl_extend_bottom = max(0.0, cyl_extend_bottom)
+    cyl_extend_top    = max(0.0, cyl_extend_top)
+
+    z_off             = (1.0 - __metric_thread_effective_ratio()) * pitch / 4
+    t_start           = z_start
+    t_length          = length
+    if bottom_relief:
+        t_start       = t_start  + (2 * z_off)
+        t_length      = t_length - (2 * z_off)
+    if top_relief:
+        t_length      = t_length - (2 * z_off)
+    outer_r           = (force_outer_radius if (force_outer_radius > 0.0) else
+                         metric_thread_major_radius(diameter,pitch))
+    inner_r           = metric_thread_minor_radius(diameter,pitch)
+    epsilon           = 0
+    inner_r_adj       = inner_r
+    inner_z_budge     = 0
+    if use_epsilon:
+        epsilon       = (z_off/3) / math.tan(__deg2rad(metric_thread_angle()))
+        inner_r_adj   = inner_r - epsilon
+        inner_z_budge = math.tan(__deg2rad(metric_thread_angle())) * epsilon
+
+    if envelope:
+        threads = cq.Workplane("XZ")
+        threads = threads.moveTo(inner_r_adj, -pitch)
+        threads = threads.lineTo(outer_r, -pitch)
+        threads = threads.lineTo(outer_r, t_length + pitch)
+        threads = threads.lineTo(inner_r_adj, t_length + pitch)
+        threads = threads.close()
+        threads = threads.revolve()
+
+    else: # Not envelope, cut the threads.
+        wire = cq.Wire.makeHelix(pitch=pitch*n_starts,
+                                 height=t_length+pitch,
+                                 radius=inner_r)
+        wire = wire.translate((0,0,-pitch/2))
+        wire = wire.rotate(startVector=(0,0,0), endVector=(0,0,1),
+                           angleDegrees=360*(-pitch/2)/(pitch*n_starts))
+        d_mid = ((metric_thread_major_radius(diameter,pitch) - outer_r)
+                 * math.tan(__deg2rad(metric_thread_angle())))
+        thread = cq.Workplane("XZ")
+        thread = thread.moveTo(inner_r_adj, -pitch/2 + z_off - inner_z_budge)
+        thread = thread.lineTo(outer_r, -(z_off + d_mid))
+        thread = thread.lineTo(outer_r, z_off + d_mid)
+        thread = thread.lineTo(inner_r_adj, pitch/2 - z_off + inner_z_budge)
+        thread = thread.close()
+        thread = thread.sweep(wire, isFrenet=True)
+        threads = thread
+        for addl_start in range(1, n_starts):
+            # TODO: Incremental/cumulative rotation may not be as accurate as
+            #       keeping 'thread' intact and rotating it by correct amount
+            #       on each iteration.  However, changing the code in that
+            #       regard may disrupt the delicate nature of workarounds
+            #       with repsect to quirks in the underlying B-rep library.
+            thread = thread.rotate(axisStartPoint=(0,0,0),
+                                   axisEndPoint=(0,0,1),
+                                   angleDegrees=360/n_starts)
+            threads = threads.union(thread)
+
+    square_shave = cq.Workplane("XY")
+    square_shave = square_shave.box(length=outer_r*3, width=outer_r*3,
+                                    height=pitch*2, centered=True)
+    square_shave = square_shave.translate((0,0,-pitch)) # Because centered.
+    # Always cut the top and bottom square.  Otherwise things don't play nice.
+    threads = threads.cut(square_shave)
+
+    if bottom_lead_in:
+        delta_r = outer_r - inner_r
+        rise = math.tan(__deg2rad(metric_thread_angle())) * delta_r
+        lead_in = cq.Workplane("XZ")
+        lead_in = lead_in.moveTo(inner_r - delta_r, -rise)
+        lead_in = lead_in.lineTo(outer_r + delta_r, 2 * rise)
+        lead_in = lead_in.lineTo(outer_r + delta_r, -pitch - rise)
+        lead_in = lead_in.lineTo(inner_r - delta_r, -pitch - rise)
+        lead_in = lead_in.close()
+        lead_in = lead_in.revolve()
+        threads = threads.cut(lead_in)
+
+    # This was originally a workaround to the anomalous B-rep computation where
+    # the top of base cylinder is flush with top of threads, without the use of
+    # lead-in.  It turns out that preferring the use of the 'render_cyl_early'
+    # strategy alleviates other problems as well.
+    render_cyl_early = (base_cylinder and ((not top_relief) and
+                                           (not (cyl_extend_top > 0.0)) and
+                                           (not envelope)))
+    render_cyl_late = (base_cylinder and (not render_cyl_early))
+    if render_cyl_early:
+        cyl = cq.Workplane("XY")
+        cyl = cyl.circle(radius=inner_r)
+        cyl = cyl.extrude(until=length+pitch+cyl_extend_bottom)
+        # Make rotation of cylinder consistent with non-workaround case.
+        cyl = cyl.rotate(axisStartPoint=(0,0,0), axisEndPoint=(0,0,1),
+                         angleDegrees=-(360*t_start/(pitch*n_starts)))
+        cyl = cyl.translate((0,0,-t_start+(z_start-cyl_extend_bottom)))
+        threads = threads.union(cyl)
+
+    # Next, make cuts at the top.
+    square_shave = square_shave.translate((0,0,pitch*2+t_length))
+    threads = threads.cut(square_shave)
+
+    if top_lead_in:
+        delta_r = outer_r - inner_r
+        rise = math.tan(__deg2rad(metric_thread_angle())) * delta_r
+        lead_in = cq.Workplane("XZ")
+        lead_in = lead_in.moveTo(inner_r - delta_r, t_length + rise)
+        lead_in = lead_in.lineTo(outer_r + delta_r, t_length - (2 * rise))
+        lead_in = lead_in.lineTo(outer_r + delta_r, t_length + pitch + rise)
+        lead_in = lead_in.lineTo(inner_r - delta_r, t_length + pitch + rise)
+        lead_in = lead_in.close()
+        lead_in = lead_in.revolve()
+        threads = threads.cut(lead_in)
+
+    # Place the threads into position.
+    threads = threads.translate((0,0,t_start))
+    if (not envelope):
+        threads = threads.rotate(axisStartPoint=(0,0,0), axisEndPoint=(0,0,1),
+                                 angleDegrees=360*t_start/(pitch*n_starts))
+
+    if render_cyl_late:
+        cyl = cq.Workplane("XY")
+        cyl = cyl.circle(radius=inner_r)
+        cyl = cyl.extrude(until=length+cyl_extend_bottom+cyl_extend_top)
+        cyl = cyl.translate((0,0,z_start-cyl_extend_bottom))
+        threads = threads.union(cyl)
+
+    return threads
+
+
+###############################################################################
+# Generates female threads without a base tube, unless 'base_tube_od' is set to
+# something which is sufficiently greater than 'diameter' parameter.  Please
+# note that 'use_epsilon' is activated by default, which causes a slight budge
+# in the major radius, outwards, so that overlaps would be created with outer
+# tubes.  (Does not affect thread profile inside of tube or beyond extents.)
+###############################################################################
+def internal_metric_thread(diameter,  # Required parameter, e.g. 3.0 for M3x0.5
+                           pitch,     # Required parameter, e.g. 0.5 for M3x0.5
+                           length,    # Required parameter, e.g. 2.0.
+                           z_start=0.0,
+                           n_starts=1,
+                           bottom_chamfer=False, # Chamfer is at same angle as
+                           top_chamfer   =False, # thread, namely 30 degrees.
+                           bottom_relief=False, # Add relief groove to start or
+                           top_relief   =False, # end of threads (shorten).
+                           use_epsilon=True, # For outer cylinder overlap.
+                           # The base tube outer diameter must be sufficiently
+                           # large for tube to be rendered.  Otherwise ignored.
+                           base_tube_od=-1.0,
+                           tube_extend_bottom=-1.0,
+                           tube_extend_top=-1.0,
+                           envelope=False): # Draw only envelope, don't cut.
+
+    tube_extend_bottom = max(0.0, tube_extend_bottom)
+    tube_extend_top    = max(0.0, tube_extend_top)
+
+    z_off              = (1.0 - __metric_thread_effective_ratio()) * pitch / 4
+    t_start            = z_start
+    t_length           = length
+    if bottom_relief:
+        t_start        = t_start  + (2 * z_off)
+        t_length       = t_length - (2 * z_off)
+    if top_relief:
+        t_length       = t_length - (2 * z_off)
+    outer_r            = metric_thread_major_radius(diameter,pitch,
+                                                    internal=True)
+    inner_r            = metric_thread_minor_radius(diameter,pitch,
+                                                    internal=True)
+    epsilon            = 0
+    outer_r_adj        = outer_r
+    outer_z_budge      = 0
+    if use_epsilon:
+        # High values of 'epsilon' sometimes cause entire starts to disappear.
+        epsilon        = (z_off/5) / math.tan(__deg2rad(metric_thread_angle()))
+        outer_r_adj    = outer_r + epsilon
+        outer_z_budge  = math.tan(__deg2rad(metric_thread_angle())) * epsilon
+
+    if envelope:
+        threads = cq.Workplane("XZ")
+        threads = threads.moveTo(outer_r_adj, -pitch)
+        threads = threads.lineTo(inner_r, -pitch)
+        threads = threads.lineTo(inner_r, t_length + pitch)
+        threads = threads.lineTo(outer_r_adj, t_length + pitch)
+        threads = threads.close()
+        threads = threads.revolve()
+
+    else: # Not envelope, cut the threads.
+        wire = cq.Wire.makeHelix(pitch=pitch*n_starts,
+                                 height=t_length+pitch,
+                                 radius=inner_r)
+        wire = wire.translate((0,0,-pitch/2))
+        wire = wire.rotate(startVector=(0,0,0), endVector=(0,0,1),
+                           angleDegrees=360*(-pitch/2)/(pitch*n_starts))
+        thread = cq.Workplane("XZ")
+        thread = thread.moveTo(outer_r_adj, -pitch/2 + z_off - outer_z_budge)
+        thread = thread.lineTo(inner_r, -z_off)
+        thread = thread.lineTo(inner_r, z_off)
+        thread = thread.lineTo(outer_r_adj, pitch/2 - z_off + outer_z_budge)
+        thread = thread.close()
+        thread = thread.sweep(wire, isFrenet=True)
+        threads = thread
+        for addl_start in range(1, n_starts):
+            # TODO: Incremental/cumulative rotation may not be as accurate as
+            #       keeping 'thread' intact and rotating it by correct amount
+            #       on each iteration.  However, changing the code in that
+            #       regard may disrupt the delicate nature of workarounds
+            #       with repsect to quirks in the underlying B-rep library.
+            thread = thread.rotate(axisStartPoint=(0,0,0),
+                                   axisEndPoint=(0,0,1),
+                                   angleDegrees=360/n_starts)
+            threads = threads.union(thread)
+        # Rotate so that the external threads would align.
+        threads = threads.rotate(axisStartPoint=(0,0,0), axisEndPoint=(0,0,1),
+                                 angleDegrees=180/n_starts)
+
+    square_len = max(outer_r*3, base_tube_od*1.125)
+    square_shave = cq.Workplane("XY")
+    square_shave = square_shave.box(length=square_len, width=square_len,
+                                    height=pitch*2, centered=True)
+    square_shave = square_shave.translate((0,0,-pitch)) # Because centered.
+    # Always cut the top and bottom square.  Otherwise things don't play nice.
+    threads = threads.cut(square_shave)
+
+    if bottom_chamfer:
+        delta_r = outer_r - inner_r
+        rise = math.tan(__deg2rad(metric_thread_angle())) * delta_r
+        chamfer = cq.Workplane("XZ")
+        chamfer = chamfer.moveTo(inner_r - delta_r, 2 * rise)
+        chamfer = chamfer.lineTo(outer_r + delta_r, -rise)
+        chamfer = chamfer.lineTo(outer_r + delta_r, -pitch - rise)
+        chamfer = chamfer.lineTo(inner_r - delta_r, -pitch - rise)
+        chamfer = chamfer.close()
+        chamfer = chamfer.revolve()
+        threads = threads.cut(chamfer)
+
+    # This was originally a workaround to the anomalous B-rep computation where
+    # the top of base tube is flush with top of threads w/o the use of chamfer.
+    # This is now being made consistent with the 'render_cyl_early' strategy in
+    # external_metric_thread() whereby we prefer the "render early" plan of
+    # action even in cases where a top chamfer or lead-in is used.
+    render_tube_early = ((base_tube_od > (outer_r * 2)) and
+                         (not top_relief) and
+                         (not (tube_extend_top > 0.0)) and
+                         (not envelope))
+    render_tube_late = ((base_tube_od > (outer_r * 2)) and
+                        (not render_tube_early))
+    if render_tube_early:
+        tube = cq.Workplane("XY")
+        tube = tube.circle(radius=base_tube_od/2)
+        tube = tube.circle(radius=outer_r)
+        tube = tube.extrude(until=length+pitch+tube_extend_bottom)
+        # Make rotation of cylinder consistent with non-workaround case.
+        tube = tube.rotate(axisStartPoint=(0,0,0), axisEndPoint=(0,0,1),
+                           angleDegrees=-(360*t_start/(pitch*n_starts)))
+        tube = tube.translate((0,0,-t_start+(z_start-tube_extend_bottom)))
+        threads = threads.union(tube)
+
+    # Next, make cuts at the top.
+    square_shave = square_shave.translate((0,0,pitch*2+t_length))
+    threads = threads.cut(square_shave)
+
+    if top_chamfer:
+        delta_r = outer_r - inner_r
+        rise = math.tan(__deg2rad(metric_thread_angle())) * delta_r
+        chamfer = cq.Workplane("XZ")
+        chamfer = chamfer.moveTo(inner_r - delta_r, t_length - (2 * rise))
+        chamfer = chamfer.lineTo(outer_r + delta_r, t_length + rise)
+        chamfer = chamfer.lineTo(outer_r + delta_r, t_length + pitch + rise)
+        chamfer = chamfer.lineTo(inner_r - delta_r, t_length + pitch + rise)
+        chamfer = chamfer.close()
+        chamfer = chamfer.revolve()
+        threads = threads.cut(chamfer)
+
+    # Place the threads into position.
+    threads = threads.translate((0,0,t_start))
+    if (not envelope):
+        threads = threads.rotate(axisStartPoint=(0,0,0), axisEndPoint=(0,0,1),
+                                 angleDegrees=360*t_start/(pitch*n_starts))
+
+    if render_tube_late:
+        tube = cq.Workplane("XY")
+        tube = tube.circle(radius=base_tube_od/2)
+        tube = tube.circle(radius=outer_r)
+        tube = tube.extrude(until=length+tube_extend_bottom+tube_extend_top)
+        tube = tube.translate((0,0,z_start-tube_extend_bottom))
+        threads = threads.union(tube)
+
+    return threads
diff --git a/nhf/parts/springs.py b/nhf/parts/springs.py
new file mode 100644
index 0000000..31fcf51
--- /dev/null
+++ b/nhf/parts/springs.py
@@ -0,0 +1,50 @@
+import math
+import cadquery as Cq
+
+def torsion_spring(radius=12,
+                   height=20,
+                   thickness=2,
+                   omega=90,
+                   tail_length=25):
+    """
+    Produces a torsion spring with abridged geometry since sweep is very slow in
+    cq-editor.
+    """
+    base = (
+        Cq.Workplane('XY')
+        .cylinder(height=height, radius=radius,
+                  centered=(True, True, False))
+    )
+    base.faces(">Z").tag("top")
+    base.faces("<Z").tag("bot")
+    result = (
+        base
+        .cylinder(height=height, radius=radius - thickness, combine='s',
+                  centered=(True, True, True))
+        .transformed(
+            offset=(0, radius-thickness),
+            rotate=(0, 0, 0))
+        .box(
+            length=tail_length,
+            width=thickness,
+            height=thickness,
+            centered=False)
+        .copyWorkplane(Cq.Workplane('XY'))
+        .transformed(
+            offset=(0, 0, height - thickness),
+            rotate=(0, 0, omega))
+        .center(-tail_length, radius-thickness)
+        .box(
+            length=tail_length,
+            width=thickness,
+            height=thickness,
+            centered=False)
+    )
+    result.polyline([(0, radius, 0), (tail_length, radius, 0)],
+                    forConstruction=True).tag("directrix_bot")
+    c, s = math.cos(omega * math.pi / 180), math.sin(omega * math.pi / 180)
+    result.polyline([
+        (s * tail_length, c * radius - s * tail_length, height),
+        (c * tail_length + s * radius, c * radius - s * tail_length, height)],
+        forConstruction=True).tag("directrix_top")
+    return result
diff --git a/nhf/parts/test.py b/nhf/parts/test.py
new file mode 100644
index 0000000..1e70ba2
--- /dev/null
+++ b/nhf/parts/test.py
@@ -0,0 +1,57 @@
+import unittest
+import cadquery as Cq
+from nhf.checks import binary_intersection
+from nhf.parts import joints, handle, metric_threads
+
+class TestJoints(unittest.TestCase):
+
+    def test_joint_hirth(self):
+        j = joints.HirthJoint()
+        obj = j.generate()
+        self.assertIsInstance(
+            obj.val().solids(), Cq.Solid,
+            msg="Hirth joint must be in one piece")
+
+    def test_joints_hirth_assembly(self):
+        for n_tooth in [16, 20, 24]:
+            with self.subTest(n_tooth=n_tooth):
+                j = joints.HirthJoint()
+                assembly = j.assembly()
+                isect = binary_intersection(assembly)
+                self.assertLess(isect.Volume(), 1e-6,
+                                "Hirth joint assembly must not have intersection")
+    def test_joints_comma_assembly(self):
+        joints.comma_assembly()
+    def test_torsion_joint(self):
+        j = joints.TorsionJoint()
+        assembly = j.rider_track_assembly()
+        bbox = assembly.toCompound().BoundingBox()
+        self.assertAlmostEqual(bbox.zlen, j.total_height)
+
+
+class TestHandle(unittest.TestCase):
+
+    def test_handle_assembly(self):
+        h = handle.Handle()
+        assembly = h.connector_insertion_assembly()
+        bbox = assembly.toCompound().BoundingBox()
+        self.assertAlmostEqual(bbox.xlen, h.diam)
+        self.assertAlmostEqual(bbox.ylen, h.diam)
+        assembly = h.connector_one_side_insertion_assembly()
+        bbox = assembly.toCompound().BoundingBox()
+        self.assertAlmostEqual(bbox.xlen, h.diam)
+        self.assertAlmostEqual(bbox.ylen, h.diam)
+
+
+class TestMetricThreads(unittest.TestCase):
+
+    def test_major_radius(self):
+        major = 3.0
+        t = metric_threads.external_metric_thread(major, 0.5, 4.0, z_start= -0.85, top_lead_in=True)
+        bbox = t.val().BoundingBox()
+        self.assertAlmostEqual(bbox.xlen, major, places=3)
+        self.assertAlmostEqual(bbox.ylen, major, places=3)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/nhf/primitive.py b/nhf/primitive.py
deleted file mode 100644
index 35e1565..0000000
--- a/nhf/primitive.py
+++ /dev/null
@@ -1,12 +0,0 @@
-import cadquery as Cq
-
-def mystery():
-    return (
-        Cq.Workplane("XY")
-        .box(10, 5, 5)
-        .faces(">Z")
-        .workplane()
-        .hole(1)
-        .edges("|Z")
-        .fillet(2)
-    )
diff --git a/nhf/test.py b/nhf/test.py
new file mode 100644
index 0000000..c6e94a6
--- /dev/null
+++ b/nhf/test.py
@@ -0,0 +1,24 @@
+import unittest
+import cadquery as Cq
+from nhf.build import Model, target
+
+class BuildScaffold(Model):
+
+    @target(name="obj1")
+    def o1(self):
+        return Cq.Solid.makeBox(10, 10, 10)
+
+    def o2(self):
+        return Cq.Solid.makeCylinder(10, 20)
+
+class TestBuild(unittest.TestCase):
+
+    def test_build_scaffold(self):
+        s = BuildScaffold()
+        names = ["obj1"]
+        self.assertEqual(s.target_names, names)
+        self.assertEqual(s.check_all(), len(names))
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/nhf/touhou/__init__.py b/nhf/touhou/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/nhf/touhou/houjuu_nue/README.org b/nhf/touhou/houjuu_nue/README.org
new file mode 100644
index 0000000..1649107
--- /dev/null
+++ b/nhf/touhou/houjuu_nue/README.org
@@ -0,0 +1,14 @@
+#+title: Cosplay: Houjuu Nue
+
+* Controller
+
+This part describes the electrical connections and the microcontroller code.
+
+* Structure
+
+This part describes the 3d printed and laser cut structures. ~structure.blend~
+is an overall sketch of the shapes and looks of the wing.
+
+* Pattern
+
+This part describes the sewing patterns.
diff --git a/nhf/touhou/houjuu_nue/__init__.py b/nhf/touhou/houjuu_nue/__init__.py
new file mode 100644
index 0000000..fe26116
--- /dev/null
+++ b/nhf/touhou/houjuu_nue/__init__.py
@@ -0,0 +1,345 @@
+"""
+To build, execute
+```
+python3 nhf/touhou/houjuu_nue/__init__.py
+```
+
+This cosplay consists of 3 components:
+
+## Trident
+
+The trident is composed of individual segments, made of acrylic, and a 3D
+printed head (convention rule prohibits metal) with a metallic paint. To ease
+transportation, the trident handle has individual segments with threads and can
+be assembled on site.
+
+## Snake
+
+A 3D printed snake with a soft material so it can wrap around and bend
+
+## Wings
+
+This is the crux of the cosplay and the most complex component. The wings mount
+on a wearable harness. Each wing consists of 4 segments with 3 joints. Parts of
+the wing which demands transluscency are created from 1/16" acrylic panels.
+These panels serve double duty as the exoskeleton.
+
+The wings are labeled r1,r2,r3,l1,l2,l3. The segments of the wings are labeled
+from root to tip s0 (root),
+s1, s2, s3. The joints are named (from root to tip)
+shoulder, elbow, wrist in analogy with human anatomy.
+"""
+from dataclasses import dataclass, field
+import unittest
+import cadquery as Cq
+from nhf import Material, Role
+from nhf.build import Model, TargetKind, target
+from nhf.parts.joints import HirthJoint
+from nhf.parts.handle import Handle
+import nhf.touhou.houjuu_nue.wing as MW
+import nhf.touhou.houjuu_nue.trident as MT
+
+@dataclass
+class Parameters(Model):
+    """
+    Defines dimensions for the Houjuu Nue cosplay
+    """
+
+    # Thickness of the exoskeleton panel in millimetres
+    panel_thickness: float = 25.4 / 16
+
+    # Harness
+    harness_thickness: float = 25.4 / 8
+    harness_width: float = 300
+    harness_height: float = 400
+    harness_fillet: float = 10
+
+    harness_wing_base_pos: list[tuple[str, float, float]] = field(default_factory=lambda: [
+        ("r1", 70, 150),
+        ("l1", -70, 150),
+        ("r2", 100, 0),
+        ("l2", -100, 0),
+        ("r3", 70, -150),
+        ("l3", -70, -150),
+    ])
+
+    # Holes drilled onto harness for attachment with HS joint
+    harness_to_root_conn_diam: float = 6
+
+    hs_hirth_joint: HirthJoint = field(default_factory=lambda: HirthJoint(
+        radius=30,
+        radius_inner=20,
+        tooth_height=10,
+        base_height=5
+    ))
+
+    # Wing root properties
+    #
+    # The Houjuu-Scarlett joint mechanism at the base of the wing
+    hs_joint_base_width: float = 85
+    hs_joint_base_thickness: float = 10
+    hs_joint_corner_fillet: float = 5
+    hs_joint_corner_cbore_diam: float = 12
+    hs_joint_corner_cbore_depth: float = 2
+    hs_joint_corner_inset: float = 12
+
+    hs_joint_axis_diam: float = 12
+    hs_joint_axis_cbore_diam: float = 20
+    hs_joint_axis_cbore_depth: float = 3
+
+    # Exterior radius of the wing root assembly
+    wing_root_radius: float = 40
+
+    """
+    Heights for various wing joints, where the numbers start from the first joint.
+    """
+    wing_r1_height: float = 100
+    wing_r1_width: float = 400
+    wing_r2_height: float = 100
+    wing_r3_height: float = 100
+
+    trident_handle: Handle = field(default_factory=lambda: Handle(
+        diam=38,
+        diam_inner=38-2 * 25.4/8,
+        # M27-3
+        diam_threading=27,
+        thread_pitch=3,
+        diam_connector_internal=18,
+        simplify_geometry=False,
+    ))
+
+    def __post_init__(self):
+        super().__init__(name="houjuu-nue")
+        assert self.wing_root_radius > self.hs_hirth_joint.radius,\
+            "Wing root must be large enough to accomodate joint"
+
+    @target(name="trident/handle-connector")
+    def handle_connector(self):
+        return self.trident_handle.connector()
+    @target(name="trident/handle-insertion")
+    def handle_insertion(self):
+        return self.trident_handle.insertion()
+
+
+    def harness_profile(self) -> Cq.Sketch:
+        """
+        Creates the harness shape
+        """
+        w, h = self.harness_width / 2, self.harness_height / 2
+        sketch = (
+            Cq.Sketch()
+            .polygon([
+                (0.7 * w, h),
+                (w, 0),
+                (0.7 * w, -h),
+                (0.7 * -w, -h),
+                (-w, 0),
+                (0.7 * -w, h),
+            ])
+            #.rect(self.harness_width, self.harness_height)
+            .vertices()
+            .fillet(self.harness_fillet)
+        )
+        for tag, x, y in self.harness_wing_base_pos:
+            conn = [(px + x, py + y) for px, py in self.hs_joint_harness_conn()]
+            sketch = (
+                sketch
+                .push(conn)
+                .tag(tag)
+                .circle(self.harness_to_root_conn_diam / 2, mode='s')
+                .reset()
+            )
+        return sketch
+
+    @target(name="harness", kind=TargetKind.DXF)
+    def harness(self) -> Cq.Shape:
+        """
+        Creates the harness shape
+        """
+        result = (
+            Cq.Workplane('XZ')
+            .placeSketch(self.harness_profile())
+            .extrude(self.harness_thickness)
+        )
+        result.faces(">Y").tag("mount")
+        plane = result.faces(">Y").workplane()
+        for tag, x, y in self.harness_wing_base_pos:
+            conn = [(px + x, py + y) for px, py in self.hs_joint_harness_conn()]
+            for i, (px, py) in enumerate(conn):
+                (
+                    plane
+                    .moveTo(px, py)
+                    .circle(1, forConstruction='True')
+                    .edges()
+                    .tag(f"{tag}_{i}")
+                )
+        return result
+
+    def hs_joint_harness_conn(self) -> list[tuple[int, int]]:
+        """
+        Generates a set of points corresponding to the connectorss
+        """
+        dx = self.hs_joint_base_width / 2 - self.hs_joint_corner_inset
+        return [
+            (dx, dx),
+            (dx, -dx),
+            (-dx, -dx),
+            (-dx, dx),
+        ]
+
+    @target(name="hs_joint_parent")
+    def hs_joint_parent(self):
+        """
+        Parent part of the Houjuu-Scarlett joint, which is composed of a Hirth
+        coupling, a cylindrical base, and a mounting base.
+        """
+        hirth = self.hs_hirth_joint.generate()
+        conn = self.hs_joint_harness_conn()
+        result = (
+            Cq.Workplane('XY')
+            .box(
+                self.hs_joint_base_width,
+                self.hs_joint_base_width,
+                self.hs_joint_base_thickness,
+                centered=(True, True, False))
+            .translate((0, 0, -self.hs_joint_base_thickness))
+            .edges("|Z")
+            .fillet(self.hs_joint_corner_fillet)
+            .faces(">Z")
+            .workplane()
+            .pushPoints(conn)
+            .cboreHole(
+                diameter=self.harness_to_root_conn_diam,
+                cboreDiameter=self.hs_joint_corner_cbore_diam,
+                cboreDepth=self.hs_joint_corner_cbore_depth)
+        )
+        # Creates a plane parallel to the holes but shifted to the base
+        plane = result.faces(">Z").workplane(offset=-self.hs_joint_base_thickness)
+
+        for i, (px, py) in enumerate(conn):
+            (
+                plane
+                .pushPoints([(px, py)])
+                .circle(1, forConstruction='True')
+                .edges()
+                .tag(f"h{i}")
+            )
+        result = (
+            result
+            .faces(">Z")
+            .workplane()
+            .union(hirth, tol=0.1)
+            .clean()
+        )
+        result = (
+            result.faces("<Z")
+            .workplane()
+            .cboreHole(
+                diameter=self.hs_joint_axis_diam,
+                cboreDiameter=self.hs_joint_axis_cbore_diam,
+                cboreDepth=self.hs_joint_axis_cbore_depth,
+            )
+            .clean()
+        )
+        result.faces("<Z").tag("base")
+        return result
+
+    @target(name="wing_root")
+    def wing_root(self) -> Cq.Assembly:
+        """
+        Generate the wing root which contains a Hirth joint at its base and a
+        rectangular opening on its side, with the necessary interfaces.
+        """
+        return MW.wing_root(joint=self.hs_hirth_joint)
+
+    def wing_r1_profile(self) -> Cq.Sketch:
+        """
+        Generates the first wing segment profile, with the wing root pointing in
+        the positive x axis.
+        """
+        # Depression of the wing middle
+        bend = 200
+        factor = 0.7
+        result = (
+            Cq.Sketch()
+            .segment((0, 0), (0, self.wing_r1_height))
+            .spline([
+                (0, self.wing_r1_height),
+                (0.5 * self.wing_r1_width, self.wing_r1_height - factor * bend),
+                (self.wing_r1_width, self.wing_r1_height - bend),
+            ])
+            .segment(
+                (self.wing_r1_width, self.wing_r1_height - bend),
+                (self.wing_r1_width, -bend),
+            )
+            .spline([
+                (self.wing_r1_width, - bend),
+                (0.5 * self.wing_r1_width, - factor * bend),
+                (0, 0),
+            ])
+            .assemble()
+        )
+        return result
+
+    def wing_r1(self) -> Cq.Solid:
+        profile = self.wing_r1_profile()
+        result = (
+            Cq.Workplane("XY")
+            .placeSketch(profile)
+            .extrude(self.panel_thickness)
+            .val()
+        )
+        return result
+
+    ######################
+    #    Assemblies      #
+    ######################
+
+    def trident_assembly(self):
+        return MT.trident_assembly(self.trident_handle)
+
+    def harness_assembly(self):
+        harness = self.harness()
+        result = (
+            Cq.Assembly()
+            .add(harness, name="base", color=Material.WOOD_BIRCH.color)
+            .constrain("base", "Fixed")
+        )
+        for name in ["l1", "l2", "l3", "r1", "r2", "r3"]:
+            j = self.hs_joint_parent()
+            (
+                result
+                .add(j, name=name, color=Role.PARENT.color)
+                .constrain("base?mount", f"{name}?base", "Axis")
+            )
+            for i in range(4):
+                result.constrain(f"base?{name}_{i}", f"{name}?h{i}", "Point")
+        result.solve()
+        return result
+
+    def wings_assembly(self):
+        """
+        Assembly of harness with all the wings
+        """
+        a_tooth = self.hs_hirth_joint.tooth_angle
+
+        result = (
+            Cq.Assembly()
+            .add(self.harness_assembly(), name="harness", loc=Cq.Location((0, 0, 0)))
+            .add(self.wing_root(), name="w0_r1")
+            .add(self.wing_root(), name="w0_l1")
+            .constrain("harness/base", "Fixed")
+            .constrain("w0_r1/joint?mate", "harness/r1?mate", "Plane")
+            .constrain("w0_r1/joint?directrix", "harness/r1?directrix",
+                       "Axis", param=7 * a_tooth)
+            .constrain("w0_l1/joint?mate", "harness/l1?mate", "Plane")
+            .constrain("w0_l1/joint?directrix", "harness/l1?directrix",
+                       "Axis", param=-1 * a_tooth)
+            .solve()
+        )
+        return result
+
+
+if __name__ == '__main__':
+    p = Parameters()
+    p.build_all()
diff --git a/nhf/touhou/houjuu_nue/controller/controller.ino b/nhf/touhou/houjuu_nue/controller/controller.ino
new file mode 100644
index 0000000..245d39d
--- /dev/null
+++ b/nhf/touhou/houjuu_nue/controller/controller.ino
@@ -0,0 +1,68 @@
+#include <FastLED.h>
+
+// Main LED strip setup
+#define LED_PIN 5
+#define NUM_LEDS 100
+#define LED_PART 50
+#define BRIGHTNESS 250
+#define LED_TYPE WS2811
+CRGB leds[NUM_LEDS];
+
+CRGB color_red;
+CRGB color_blue;
+CRGB color_green;
+
+#define DIAG_PIN 6
+
+
+void setup() {
+  // Calculate colors
+  hsv2rgb_spectrum(CHSV(4, 255, 100), color_red);
+  hsv2rgb_spectrum(CHSV(170, 255, 100), color_blue);
+  hsv2rgb_spectrum(CHSV(90, 255, 100), color_green);
+  pinMode(LED_BUILTIN, OUTPUT);
+  pinMode(LED_PIN, OUTPUT);
+  pinMode(DIAG_PIN, OUTPUT);
+  
+  // Main LED strip
+  FastLED.addLeds<LED_TYPE, LED_PIN, RGB>(leds, NUM_LEDS);
+}
+
+void loop() {
+  fill_segmented(CRGB::Green, CRGB::Orange);
+  delay(500);
+
+  flash(leds, NUM_LEDS, color_red, 10, 20);
+  delay(500);
+  flash(leds, NUM_LEDS, color_blue, 10, 20);
+  delay(500);
+}
+
+void fill_segmented(CRGB c1, CRGB c2)
+{
+  //fill_solid(leds, LED_PART, c1);
+  fill_gradient_RGB(leds, LED_PART, CRGB::Black ,c1);
+  fill_gradient_RGB(leds + LED_PART, NUM_LEDS - LED_PART, CRGB::Black, c2);
+  FastLED.show(); 
+}
+void flash(CRGB *ptr, uint16_t num, CRGB const& lead, int steps, int step_time)
+{
+  digitalWrite(LED_BUILTIN, LOW);
+
+  //fill_solid(leds, NUM_LEDS, CRGB::Black);
+  for (int i = 0; i < steps; ++i)
+  {
+    uint8_t factor = 255 * i / steps;
+    analogWrite(DIAG_PIN, factor);
+    CRGB tail = blend(lead, CRGB::Black, factor);
+    uint16_t front = factor * (int) num / 255;
+    fill_solid(ptr, front, tail);
+    //fill_gradient_RGB(ptr, front, tail, lead);
+    //fill_solid(leds + front, NUM_LEDS - front, CRGB::Black);
+    FastLED.show(); 
+    delay(step_time);
+  }
+  fill_gradient_RGB(ptr, num, CRGB::Black, lead);
+  FastLED.show(); 
+  analogWrite(DIAG_PIN, LOW);
+}
\ No newline at end of file
diff --git a/nhf/touhou/houjuu_nue/test.py b/nhf/touhou/houjuu_nue/test.py
new file mode 100644
index 0000000..0b63ecd
--- /dev/null
+++ b/nhf/touhou/houjuu_nue/test.py
@@ -0,0 +1,39 @@
+import unittest
+import cadquery as Cq
+import nhf.touhou.houjuu_nue as M
+
+class Test(unittest.TestCase):
+
+    def test_hs_joint_parent(self):
+        p = M.Parameters()
+        obj = p.hs_joint_parent()
+        self.assertIsInstance(obj.val().solids(), Cq.Solid, msg="H-S joint must be in one piece")
+    def test_wing_root(self):
+        p = M.Parameters()
+        obj = p.wing_root()
+        #self.assertIsInstance(obj.solids(), Cq.Solid, msg="Wing root must be in one piece")
+        bbox = obj.val().BoundingBox()
+
+        msg = "Must fix 256^3 bbox"
+        self.assertLess(bbox.xlen, 255, msg=msg)
+        self.assertLess(bbox.ylen, 255, msg=msg)
+        self.assertLess(bbox.zlen, 255, msg=msg)
+    def test_wing_root(self):
+        p = M.Parameters()
+        p.wing_root()
+    def test_wings_assembly(self):
+        p = M.Parameters()
+        p.wings_assembly()
+    def test_harness_assembly(self):
+        p = M.Parameters()
+        p.harness_assembly()
+    def test_trident_assembly(self):
+        p = M.Parameters()
+        assembly = p.trident_assembly()
+        bbox = assembly.toCompound().BoundingBox()
+        length = bbox.zlen
+        self.assertGreater(length, 1300)
+        self.assertLess(length, 1700)
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/nhf/touhou/houjuu_nue/trident.py b/nhf/touhou/houjuu_nue/trident.py
new file mode 100644
index 0000000..312ff85
--- /dev/null
+++ b/nhf/touhou/houjuu_nue/trident.py
@@ -0,0 +1,42 @@
+import math
+import cadquery as Cq
+from nhf import Material
+from nhf.parts.handle import Handle
+
+def trident_assembly(
+        handle: Handle,
+        handle_segment_length: float = 24*25.4,
+        terminal_height=100):
+    def segment():
+        return handle.segment(handle_segment_length)
+
+    terminal = (
+        handle
+        .one_side_connector(height=terminal_height)
+        .faces(">Z")
+        .hole(15, terminal_height + handle.insertion_length - 10)
+    )
+    mat_i = Material.PLASTIC_PLA
+    mat_s = Material.ACRYLIC_BLACK
+    assembly = (
+        Cq.Assembly()
+        .add(handle.insertion(), name="i0", color=mat_i.color)
+        .constrain("i0", "Fixed")
+        .add(segment(), name="s1", color=mat_s.color)
+        .constrain("i0?rim", "s1?mate1", "Plane", param=0)
+        .add(handle.insertion(), name="i1", color=mat_i.color)
+        .add(handle.connector(), name="c1", color=mat_i.color)
+        .add(handle.insertion(), name="i2", color=mat_i.color)
+        .constrain("s1?mate2", "i1?rim", "Plane", param=0)
+        .constrain("i1?mate", "c1?mate1", "Plane")
+        .constrain("i2?mate", "c1?mate2", "Plane")
+        .add(segment(), name="s2", color=mat_s.color)
+        .constrain("i2?rim", "s2?mate1", "Plane", param=0)
+        .add(handle.insertion(), name="i3", color=mat_i.color)
+        .constrain("s2?mate2", "i3?rim", "Plane", param=0)
+        .add(handle.one_side_connector(), name="head", color=mat_i.color)
+        .constrain("i3?mate", "head?mate", "Plane")
+        .add(terminal, name="terminal", color=mat_i.color)
+        .constrain("i0?mate", "terminal?mate", "Plane")
+    )
+    return assembly.solve()
diff --git a/nhf/touhou/houjuu_nue/wing.py b/nhf/touhou/houjuu_nue/wing.py
new file mode 100644
index 0000000..ea3a5c3
--- /dev/null
+++ b/nhf/touhou/houjuu_nue/wing.py
@@ -0,0 +1,226 @@
+"""
+This file describes the shapes of the wing shells. The joints are defined in
+`__init__.py`.
+"""
+import math
+import cadquery as Cq
+from nhf import Material, Role
+from nhf.parts.joints import HirthJoint
+
+def wing_root_profiles(
+        base_sweep=150,
+        wall_thickness=8,
+        base_radius=40,
+        middle_offset=30,
+        middle_height=80,
+        conn_width=40,
+        conn_height=100) -> tuple[Cq.Wire, Cq.Wire]:
+    assert base_sweep < 180
+    assert middle_offset > 0
+    theta = math.pi * base_sweep / 180
+    c, s = math.cos(theta), math.sin(theta)
+    c_1, s_1 = math.cos(theta * 0.75), math.sin(theta * 0.75)
+    c_2, s_2 = math.cos(theta / 2), math.sin(theta / 2)
+    r1 = base_radius
+    r2 = base_radius - wall_thickness
+    base = (
+        Cq.Sketch()
+        .arc(
+            (c * r1, s * r1),
+            (c_1 * r1, s_1 * r1),
+            (c_2 * r1, s_2 * r1),
+        )
+        .arc(
+            (c_2 * r1, s_2 * r1),
+            (r1, 0),
+            (c_2 * r1, -s_2 * r1),
+        )
+        .arc(
+            (c_2 * r1, -s_2 * r1),
+            (c_1 * r1, -s_1 * r1),
+            (c * r1, -s * r1),
+        )
+        .segment(
+            (c * r1, -s * r1),
+            (c * r2, -s * r2),
+        )
+        .arc(
+            (c * r2, -s * r2),
+            (c_1 * r2, -s_1 * r2),
+            (c_2 * r2, -s_2 * r2),
+        )
+        .arc(
+            (c_2 * r2, -s_2 * r2),
+            (r2, 0),
+            (c_2 * r2, s_2 * r2),
+        )
+        .arc(
+            (c_2 * r2, s_2 * r2),
+            (c_1 * r2, s_1 * r2),
+            (c * r2, s * r2),
+        )
+        .segment(
+            (c * r2, s * r2),
+            (c * r1, s * r1),
+        )
+        .assemble(tag="wire")
+        .wires().val()
+    )
+    assert isinstance(base, Cq.Wire)
+
+    # The interior sweep is given by theta, but the exterior sweep exceeds the
+    # interior sweep so the wall does not become thinner towards the edges.
+    # If the exterior sweep is theta', it has to satisfy
+    #
+    # sin(theta) * r2 + wall_thickness = sin(theta') * r1
+    x, y = conn_width / 2, middle_height / 2
+    t = wall_thickness
+    dx = middle_offset
+    middle = (
+        Cq.Sketch()
+        # Interior arc, top point
+        .arc(
+            (x - t, y - t),
+            (x - t + dx, 0),
+            (x - t, -y + t),
+        )
+        .segment(
+            (x - t, -y + t),
+            (-x, -y+t)
+        )
+        .segment((-x, -y))
+        .segment((x, -y))
+        # Outer arc, bottom point
+        .arc(
+            (x, -y),
+            (x + dx, 0),
+            (x, y),
+        )
+        .segment(
+            (x, y),
+            (-x, y)
+        )
+        .segment((-x, y-t))
+        #.segment((x2, a))
+        .close()
+        .assemble(tag="wire")
+        .wires().val()
+    )
+    assert isinstance(middle, Cq.Wire)
+
+    x, y = conn_width / 2, conn_height / 2
+    t = wall_thickness
+    tip = (
+        Cq.Sketch()
+        .segment((-x, y), (x, y))
+        .segment((x, -y))
+        .segment((-x, -y))
+        .segment((-x, -y+t))
+        .segment((x-t, -y+t))
+        .segment((x-t, y-t))
+        .segment((-x, y-t))
+        .close()
+        .assemble(tag="wire")
+        .wires().val()
+    )
+    return base, middle, tip
+
+
+def wing_root(joint: HirthJoint,
+              bolt_diam: int = 12,
+              union_tol=1e-4,
+              attach_diam=8,
+              conn_width=40,
+              conn_height=100,
+              wall_thickness=8) -> Cq.Assembly:
+    """
+    Generate the contiguous components of the root wing segment
+    """
+    tip_centre = Cq.Vector((-150, 0, -80))
+    attach_points = [
+        (15, 0),
+        (40, 0),
+    ]
+    root_profile, middle_profile, tip_profile = wing_root_profiles(
+        conn_width=conn_width,
+        conn_height=conn_height,
+        wall_thickness=8,
+    )
+    middle_profile = middle_profile.located(Cq.Location(
+        (-40, 0, -40), (0, 30, 0)
+    ))
+    antetip_profile = tip_profile.located(Cq.Location(
+        (-95, 0, -75), (0, 60, 0)
+    ))
+    tip_profile = tip_profile.located(Cq.Location(
+        tip_centre, (0, 90, 0)
+    ))
+    profiles = [
+        root_profile,
+        middle_profile,
+        antetip_profile,
+        tip_profile,
+    ]
+    result = None
+    for p1, p2 in zip(profiles[:-1], profiles[1:]):
+        seg = (
+            Cq.Workplane('XY')
+            .add(p1)
+            .toPending()
+            .workplane() # This call is necessary
+            .add(p2)
+            .toPending()
+            .loft()
+        )
+        if result:
+            result = result.union(seg, tol=union_tol)
+        else:
+            result = seg
+    result = (
+        result
+        # Create connector holes
+        .copyWorkplane(
+            Cq.Workplane('bottom', origin=tip_centre + Cq.Vector((0, -50, 0)))
+        )
+        .pushPoints(attach_points)
+        .hole(attach_diam)
+    )
+    # Generate attach point tags
+
+    for sign in [False, True]:
+        y = conn_height / 2 - wall_thickness
+        side = "bottom" if sign else "top"
+        y = y if sign else -y
+        plane = (
+            result
+            # Create connector holes
+            .copyWorkplane(
+                Cq.Workplane(side, origin=tip_centre +
+                             Cq.Vector((0, y, 0)))
+            )
+        )
+        for i, (px, py) in enumerate(attach_points):
+            (
+                plane
+                .moveTo(px, py)
+                .eachpoint(Cq.Vertex.makeVertex(0, 0, 0))
+                .tag(f"conn_{side}{i}")
+            )
+
+    result.faces("<Z").tag("base")
+    result.faces(">X").tag("conn")
+
+    j = (
+        joint.generate(is_mated=True)
+        .faces("<Z")
+        .hole(bolt_diam)
+    )
+
+    color = Material.PLASTIC_PLA.color
+    result = (
+        Cq.Assembly()
+        .add(result, name="scaffold", color=color)
+        .add(j, name="joint", color=Role.CHILD.color,
+             loc=Cq.Location((0, 0, -joint.total_height)))
+    )
+    return result
diff --git a/poetry.lock b/poetry.lock
index 69302dd..756343c 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -857,4 +857,4 @@ files = [
 [metadata]
 lock-version = "2.0"
 python-versions = "^3.10"
-content-hash = "caf46b526858dbf2960b0204782140ad072d96f7b3f161ac7e9db0d9b709b25a"
+content-hash = "ec47ccffd60fbda610a5c3725fc064a08b1b794f23084672bd62beb20b1b19f7"
diff --git a/pyproject.toml b/pyproject.toml
index 3431c0e..8ca0117 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -10,6 +10,7 @@ python = "^3.10"
 cadquery = "^2.4.0"
 build123d = "^0.5.0"
 numpy = "^1.26.4"
+colorama = "^0.4.6"
 
 [build-system]
 requires = ["poetry-core"]