#!/usr/bin/env python3
"""Two printable solids that seat together: a hollow cylinder and a disc-cone.

The cylinder is a 50 mm tube. Its outer diameter is 72 mm and the wall is
1.5 mm, so the bore is 69 mm. The bottom edge is broken by evenly spaced
half-arch cutouts. A requested 2 mm gap between those arches does not divide
the circumference into a whole number of arches, so the count is whichever
leaves a mid-wall gap closest to 2 mm.

The second solid is a 100 mm disc, 2 mm thick, with an exponential cone on
top. The cone narrows from the disc diameter to 5 mm over 45 mm, then a 5 mm
cone spike tapers to a point. Where that cone is just inside the cylinder
bore, 1.5 mm deep cutouts accept the solid feet left between the arches.
Each mating face of a cutout is 0.1 mm outside the foot.

Run this file to write both STL files next to the script.
"""

from __future__ import annotations

import math
from pathlib import Path

import cadquery as cq
from OCP.BRepClass3d import BRepClass3d_SolidClassifier
from OCP.TopAbs import TopAbs_IN
from OCP.gp import gp_Pnt

# --- Hollow cylinder dimensions, millimetres ---

# Axial length of the tube. The bottom face is z = 0 and the top face is this height.
CYLINDER_HEIGHT_MILLIMETERS = 50.0
# Driving diameter. The bore is this minus twice the wall.
CYLINDER_OUTER_DIAMETER_MILLIMETERS = 72.0
CYLINDER_WALL_THICKNESS_MILLIMETERS = 1.5
CYLINDER_INNER_DIAMETER_MILLIMETERS = (
    CYLINDER_OUTER_DIAMETER_MILLIMETERS - 2.0 * CYLINDER_WALL_THICKNESS_MILLIMETERS
)
# Radius of each half-arch cut up from the bottom edge. The opening width is twice this.
ARCH_CUTOUT_RADIUS_MILLIMETERS = 15.0
# Target clear arc, measured at mid-wall, between neighbouring arch openings.
REQUESTED_ARCH_GAP_MILLIMETERS = 2.0

CYLINDER_INNER_RADIUS_MILLIMETERS = CYLINDER_INNER_DIAMETER_MILLIMETERS / 2.0
CYLINDER_OUTER_RADIUS_MILLIMETERS = CYLINDER_OUTER_DIAMETER_MILLIMETERS / 2.0
CYLINDER_MIDWALL_RADIUS_MILLIMETERS = (
    CYLINDER_INNER_RADIUS_MILLIMETERS + CYLINDER_OUTER_RADIUS_MILLIMETERS
) / 2.0

# How far the arch cutter extends past each face of the wall, so the boolean
# goes all the way through instead of leaving a skin.
ARCH_CUTTER_AXIAL_OVERTRAVEL_MILLIMETERS = 6.0

# --- Disc, exponential cone, and spike ---

DISC_OUTER_DIAMETER_MILLIMETERS = 100.0
DISC_THICKNESS_MILLIMETERS = 2.0
# Height of the narrowing cone plus the spike, not including the disc.
CONE_AND_SPIKE_STACK_HEIGHT_MILLIMETERS = 50.0
SPIKE_BASE_DIAMETER_MILLIMETERS = 5.0
SPIKE_HEIGHT_MILLIMETERS = 5.0
# The narrowing section is the stack minus the spike: 50 mm - 5 mm.
NARROWING_CONE_HEIGHT_MILLIMETERS = (
    CONE_AND_SPIKE_STACK_HEIGHT_MILLIMETERS - SPIKE_HEIGHT_MILLIMETERS
)
# The cone ends at the spike's base diameter, then the spike continues to a point.
NARROWING_CONE_END_DIAMETER_MILLIMETERS = SPIKE_BASE_DIAMETER_MILLIMETERS

DISC_OUTER_RADIUS_MILLIMETERS = DISC_OUTER_DIAMETER_MILLIMETERS / 2.0
NARROWING_CONE_END_RADIUS_MILLIMETERS = NARROWING_CONE_END_DIAMETER_MILLIMETERS / 2.0
DISC_CONE_OVERALL_HEIGHT_MILLIMETERS = (
    DISC_THICKNESS_MILLIMETERS + NARROWING_CONE_HEIGHT_MILLIMETERS + SPIKE_HEIGHT_MILLIMETERS
)

# Radius along the cone is DISC_OUTER_RADIUS * exp(decay * height / cone_height).
# decay is negative, so the radius shrinks from the disc to the 5 mm end.
EXPONENTIAL_CONE_DECAY_CONSTANT = math.log(
    NARROWING_CONE_END_RADIUS_MILLIMETERS / DISC_OUTER_RADIUS_MILLIMETERS
)
# Straight segments used to revolve the exponential. The solid's surface is the
# chord between these samples, which sits slightly outside the true curve.
EXPONENTIAL_CONE_PROFILE_SEGMENT_COUNT = 96

# How far each arch-gap foot sinks into the cone when the cylinder is seated.
JOIN_CUTOUT_DEPTH_MILLIMETERS = 1.5
# Extra space on every mating face of a cutout so the foot can enter.
JOIN_FIT_CLEARANCE_MILLIMETERS = 0.1

# --- Output files ---

HOLLOW_CYLINDER_STL_FILENAME = "hollow_cylinder.stl"
DISC_CONE_STL_FILENAME = "disc_cone.stl"
HOLLOW_CYLINDER_STL_PATH = Path(__file__).with_name(HOLLOW_CYLINDER_STL_FILENAME)
DISC_CONE_STL_PATH = Path(__file__).with_name(DISC_CONE_STL_FILENAME)

# Mesh tolerances passed to the STL exporter.
STL_LINEAR_DEFLECTION_MILLIMETERS = 0.02
STL_ANGULAR_DEFLECTION_RADIANS = 0.05

# Arch-count search. Gaps smaller than the minimum are treated as closed.
MINIMUM_ARCH_COUNT_TO_TRY = 3
MAXIMUM_ARCH_COUNT_TO_TRY = 59
MINIMUM_OPEN_ARCH_GAP_MILLIMETERS = 0.2

# Checks that the generated solids match the dimensions above.
BOUNDING_BOX_TOLERANCE_MILLIMETERS = 0.05
POINT_CLASSIFIER_TOLERANCE_MILLIMETERS = 1e-7
# Touching seat faces can leave a hair of intersection volume. Anything larger is a clash.
MAXIMUM_SEATED_OVERLAP_VOLUME_CUBIC_MILLIMETERS = 0.5


def clear_arch_gap_arc_length_millimeters(radius_millimeters: float, arch_count: int) -> float:
    """Uncut arc between neighbouring arch openings, measured at radius_millimeters.

    Each arch is a radial cutter, so the opening's half-angle at a given
    radius is asin(arch_radius / radius). The gap is whatever is left of the
    equal angular pitch after that opening is removed.
    """
    angular_pitch_radians = 2.0 * math.pi / arch_count
    arch_opening_radians = 2.0 * math.asin(ARCH_CUTOUT_RADIUS_MILLIMETERS / radius_millimeters)
    return radius_millimeters * (angular_pitch_radians - arch_opening_radians)


def choose_arch_count() -> int:
    """Pick the arch count whose mid-wall gap is closest to the requested 2 mm.

    The circumference is not an integer multiple of (opening + requested gap),
    so the arches are spaced evenly and the gap comes out a little off 2 mm.
    """
    closest_gap_error_and_count: tuple[float, int] | None = None
    for candidate_arch_count in range(MINIMUM_ARCH_COUNT_TO_TRY, MAXIMUM_ARCH_COUNT_TO_TRY + 1):
        midwall_gap_millimeters = clear_arch_gap_arc_length_millimeters(
            CYLINDER_MIDWALL_RADIUS_MILLIMETERS,
            candidate_arch_count,
        )
        if midwall_gap_millimeters <= MINIMUM_OPEN_ARCH_GAP_MILLIMETERS:
            continue
        gap_error_millimeters = abs(midwall_gap_millimeters - REQUESTED_ARCH_GAP_MILLIMETERS)
        if (
            closest_gap_error_and_count is None
            or gap_error_millimeters < closest_gap_error_and_count[0]
        ):
            closest_gap_error_and_count = (gap_error_millimeters, candidate_arch_count)
    if closest_gap_error_and_count is None:
        raise RuntimeError("no arch count leaves a positive gap")
    return closest_gap_error_and_count[1]


def make_hollow_cylinder_wall() -> cq.Workplane:
    """Open tube before the arch cutouts. Outer circle first, then the bore."""
    return (
        cq.Workplane("XY")
        .circle(CYLINDER_OUTER_RADIUS_MILLIMETERS)
        .circle(CYLINDER_INNER_RADIUS_MILLIMETERS)
        .extrude(CYLINDER_HEIGHT_MILLIMETERS)
    )


def make_radial_arch_cutter(
    arch_index: int,
    arch_count: int,
    cutter_radius_millimeters: float,
    cutter_axis_height_millimeters: float,
    radial_start_millimeters: float,
    radial_end_millimeters: float,
) -> cq.Solid:
    """Horizontal cylinder lying on a radius, used as a half-arch cutter.

    The circle is drawn in the YZ plane and extruded along X, then rotated
    about Z to this arch's station. The axis sits on cutter_axis_height, so
    only the upper half meets a solid that begins at that height. That upper
    half is the half-arch.
    """
    arch_station_angle_degrees = 360.0 * arch_index / arch_count
    cutter_length_millimeters = radial_end_millimeters - radial_start_millimeters
    cutter_center_radius_millimeters = (radial_start_millimeters + radial_end_millimeters) / 2.0
    radial_cutter = (
        cq.Workplane("YZ")
        .circle(cutter_radius_millimeters)
        .extrude(cutter_length_millimeters / 2.0, both=True)
        .translate((cutter_center_radius_millimeters, 0.0, cutter_axis_height_millimeters))
        .rotate((0, 0, 0), (0, 0, 1), arch_station_angle_degrees)
    )
    return radial_cutter.val()


def make_cylinder_arch_cutter(arch_index: int, arch_count: int) -> cq.Solid:
    """Arch cutter for the cylinder wall, centered on the bottom plane (z = 0)."""
    radial_overtravel_millimeters = (
        CYLINDER_WALL_THICKNESS_MILLIMETERS + ARCH_CUTTER_AXIAL_OVERTRAVEL_MILLIMETERS
    ) / 2.0
    return make_radial_arch_cutter(
        arch_index,
        arch_count,
        ARCH_CUTOUT_RADIUS_MILLIMETERS,
        0.0,
        CYLINDER_MIDWALL_RADIUS_MILLIMETERS - radial_overtravel_millimeters,
        CYLINDER_MIDWALL_RADIUS_MILLIMETERS + radial_overtravel_millimeters,
    )


def build_hollow_cylinder() -> tuple[cq.Workplane, int]:
    """Wall with every half-arch removed. Returns the solid and the arch count."""
    arch_count = choose_arch_count()
    hollow_cylinder_wall = make_hollow_cylinder_wall()
    arch_cutters = cq.Compound.makeCompound(
        [make_cylinder_arch_cutter(arch_index, arch_count) for arch_index in range(arch_count)]
    )
    return hollow_cylinder_wall.cut(arch_cutters), arch_count


def point_is_inside_solid(
    solid_shape: cq.Shape,
    x_millimeters: float,
    y_millimeters: float,
    z_millimeters: float,
) -> bool:
    """True when the sample point is strictly inside the solid, not on its skin."""
    solid_classifier = BRepClass3d_SolidClassifier(
        solid_shape.wrapped,
        gp_Pnt(x_millimeters, y_millimeters, z_millimeters),
        POINT_CLASSIFIER_TOLERANCE_MILLIMETERS,
    )
    return solid_classifier.State() == TopAbs_IN


def verify_hollow_cylinder(hollow_cylinder: cq.Workplane, arch_count: int) -> None:
    """Check the tube bounds, the bore, and that every arch actually cuts the wall."""
    hollow_cylinder_shape = hollow_cylinder.val()
    if not hollow_cylinder_shape.isValid():
        raise SystemExit("hollow cylinder solid is not valid")

    bounding_box = hollow_cylinder_shape.BoundingBox()
    expected_axis_limits_millimeters = {
        "x": (bounding_box.xmin, bounding_box.xmax, -CYLINDER_OUTER_RADIUS_MILLIMETERS, CYLINDER_OUTER_RADIUS_MILLIMETERS),
        "y": (bounding_box.ymin, bounding_box.ymax, -CYLINDER_OUTER_RADIUS_MILLIMETERS, CYLINDER_OUTER_RADIUS_MILLIMETERS),
        "z": (bounding_box.zmin, bounding_box.zmax, 0.0, CYLINDER_HEIGHT_MILLIMETERS),
    }
    for axis_name, (
        actual_low_millimeters,
        actual_high_millimeters,
        expected_low_millimeters,
        expected_high_millimeters,
    ) in expected_axis_limits_millimeters.items():
        low_error_millimeters = abs(actual_low_millimeters - expected_low_millimeters)
        high_error_millimeters = abs(actual_high_millimeters - expected_high_millimeters)
        if (
            low_error_millimeters > BOUNDING_BOX_TOLERANCE_MILLIMETERS
            or high_error_millimeters > BOUNDING_BOX_TOLERANCE_MILLIMETERS
        ):
            raise SystemExit(
                f"hollow cylinder bounding box {axis_name} is "
                f"{actual_low_millimeters:.3f}..{actual_high_millimeters:.3f}"
            )

    # (x, y, z, should_be_inside, what_this_sample_is_checking)
    inside_outside_samples = [
        (CYLINDER_MIDWALL_RADIUS_MILLIMETERS, 0.0, CYLINDER_HEIGHT_MILLIMETERS / 2.0, True, "wall at mid height"),
        (CYLINDER_INNER_RADIUS_MILLIMETERS - 1.0, 0.0, CYLINDER_HEIGHT_MILLIMETERS / 2.0, False, "bore"),
        (CYLINDER_OUTER_RADIUS_MILLIMETERS + 1.0, 0.0, CYLINDER_HEIGHT_MILLIMETERS / 2.0, False, "outside"),
        (CYLINDER_MIDWALL_RADIUS_MILLIMETERS, 0.0, 2.0, False, "center of an arch"),
        (
            CYLINDER_MIDWALL_RADIUS_MILLIMETERS,
            0.0,
            ARCH_CUTOUT_RADIUS_MILLIMETERS + 2.0,
            True,
            "wall above an arch",
        ),
    ]
    # Arch 0 is on +X. The first gap center is halfway to arch 1.
    first_arch_gap_angle_radians = math.pi / arch_count
    inside_outside_samples.append(
        (
            CYLINDER_MIDWALL_RADIUS_MILLIMETERS * math.cos(first_arch_gap_angle_radians),
            CYLINDER_MIDWALL_RADIUS_MILLIMETERS * math.sin(first_arch_gap_angle_radians),
            0.4,
            True,
            "bottom rim between arches",
        )
    )
    for (
        sample_x_millimeters,
        sample_y_millimeters,
        sample_z_millimeters,
        sample_should_be_inside,
        sample_label,
    ) in inside_outside_samples:
        sample_is_inside = point_is_inside_solid(
            hollow_cylinder_shape,
            sample_x_millimeters,
            sample_y_millimeters,
            sample_z_millimeters,
        )
        if sample_is_inside != sample_should_be_inside:
            raise SystemExit(
                f"{sample_label}: expected inside={sample_should_be_inside}, got {sample_is_inside}"
            )

    arches_cut_through_wall_count = 0
    for arch_index in range(arch_count):
        arch_station_angle_radians = 2.0 * math.pi * arch_index / arch_count
        arch_center_x_millimeters = CYLINDER_MIDWALL_RADIUS_MILLIMETERS * math.cos(
            arch_station_angle_radians
        )
        arch_center_y_millimeters = CYLINDER_MIDWALL_RADIUS_MILLIMETERS * math.sin(
            arch_station_angle_radians
        )
        # Two millimetres up is still inside the half-arch, so the wall must be gone.
        if point_is_inside_solid(
            hollow_cylinder_shape,
            arch_center_x_millimeters,
            arch_center_y_millimeters,
            2.0,
        ):
            raise SystemExit(f"arch {arch_index} did not cut through the wall")
        arches_cut_through_wall_count += 1
    if arches_cut_through_wall_count != arch_count:
        raise SystemExit(f"expected {arch_count} arches, checked {arches_cut_through_wall_count}")


def exponential_cone_radius_millimeters(height_fraction: float) -> float:
    """Outer radius of the true exponential. 0 is the disc, 1 is the 5 mm end."""
    return DISC_OUTER_RADIUS_MILLIMETERS * math.exp(
        EXPONENTIAL_CONE_DECAY_CONSTANT * height_fraction
    )


def narrowing_cone_height_at_radius_millimeters(target_radius_millimeters: float) -> float:
    """Height above the disc where the revolved cone surface equals target_radius.

    The solid is a polyline, so between samples the surface is a straight chord.
    That chord sits outside the exponential. The seat is measured on the chord,
    which is the surface the cylinder actually touches.
    """
    for segment_index in range(EXPONENTIAL_CONE_PROFILE_SEGMENT_COUNT):
        segment_start_height_millimeters = (
            NARROWING_CONE_HEIGHT_MILLIMETERS
            * segment_index
            / EXPONENTIAL_CONE_PROFILE_SEGMENT_COUNT
        )
        segment_end_height_millimeters = (
            NARROWING_CONE_HEIGHT_MILLIMETERS
            * (segment_index + 1)
            / EXPONENTIAL_CONE_PROFILE_SEGMENT_COUNT
        )
        segment_start_radius_millimeters = exponential_cone_radius_millimeters(
            segment_index / EXPONENTIAL_CONE_PROFILE_SEGMENT_COUNT
        )
        segment_end_radius_millimeters = exponential_cone_radius_millimeters(
            (segment_index + 1) / EXPONENTIAL_CONE_PROFILE_SEGMENT_COUNT
        )
        radius_is_on_this_segment = (
            segment_end_radius_millimeters
            <= target_radius_millimeters
            <= segment_start_radius_millimeters
        )
        if not radius_is_on_this_segment:
            continue
        radius_span_millimeters = (
            segment_end_radius_millimeters - segment_start_radius_millimeters
        )
        if radius_span_millimeters == 0.0:
            return segment_start_height_millimeters
        segment_fraction = (
            target_radius_millimeters - segment_start_radius_millimeters
        ) / radius_span_millimeters
        return segment_start_height_millimeters + segment_fraction * (
            segment_end_height_millimeters - segment_start_height_millimeters
        )
    raise RuntimeError(f"cone never reaches radius {target_radius_millimeters:.3f} mm")


def cylinder_seat_height_millimeters() -> float:
    """World Z where the cone surface is one clearance inside the cylinder bore.

    Above this height the cone is narrower than the bore, so the tube can
    slide down over it. The arch-gap feet stop 1.5 mm below this height.
    """
    bore_radius_with_clearance_millimeters = (
        CYLINDER_INNER_RADIUS_MILLIMETERS - JOIN_FIT_CLEARANCE_MILLIMETERS
    )
    return DISC_THICKNESS_MILLIMETERS + narrowing_cone_height_at_radius_millimeters(
        bore_radius_with_clearance_millimeters
    )


def seated_cylinder_bottom_height_millimeters() -> float:
    """World Z of the cylinder's bottom face when the feet are fully in the cutouts."""
    return cylinder_seat_height_millimeters() - JOIN_CUTOUT_DEPTH_MILLIMETERS


def make_arch_gap_join_cutouts(arch_count: int) -> cq.Workplane:
    """One 1.5 mm pocket per arch gap, sized so the cylinder foot can drop in.

    The pocket is the wall annulus, grown by the fit clearance, with the arch
    openings put back using a slightly smaller arch radius. A smaller arch
    leaves a wider foot-shaped pocket, which is the 0.1 mm side clearance.
    The pocket's bottom is the floor the foot lands on.
    """
    cutout_arch_radius_millimeters = (
        ARCH_CUTOUT_RADIUS_MILLIMETERS - JOIN_FIT_CLEARANCE_MILLIMETERS
    )
    cutout_inner_radius_millimeters = (
        CYLINDER_INNER_RADIUS_MILLIMETERS - JOIN_FIT_CLEARANCE_MILLIMETERS
    )
    cutout_outer_radius_millimeters = (
        CYLINDER_OUTER_RADIUS_MILLIMETERS + JOIN_FIT_CLEARANCE_MILLIMETERS
    )
    cutout_floor_height_millimeters = seated_cylinder_bottom_height_millimeters()
    # Annulus standing on the cutout floor. The arch cutters then remove the
    # openings, leaving only the foot-shaped pockets.
    cutout_annulus = (
        cq.Workplane("XY")
        .workplane(offset=cutout_floor_height_millimeters)
        .circle(cutout_outer_radius_millimeters)
        .circle(cutout_inner_radius_millimeters)
        .extrude(JOIN_CUTOUT_DEPTH_MILLIMETERS)
    )
    radial_overtravel_millimeters = (
        CYLINDER_WALL_THICKNESS_MILLIMETERS + ARCH_CUTTER_AXIAL_OVERTRAVEL_MILLIMETERS
    ) / 2.0
    smaller_arch_cutters = [
        make_radial_arch_cutter(
            arch_index,
            arch_count,
            cutout_arch_radius_millimeters,
            cutout_floor_height_millimeters,
            CYLINDER_MIDWALL_RADIUS_MILLIMETERS - radial_overtravel_millimeters,
            CYLINDER_MIDWALL_RADIUS_MILLIMETERS + radial_overtravel_millimeters,
        )
        for arch_index in range(arch_count)
    ]
    return cutout_annulus.cut(cq.Compound.makeCompound(smaller_arch_cutters))


def make_disc_cone(arch_count: int) -> cq.Workplane:
    """Disc, exponential cone, point spike, then the arch-gap join cutouts.

    The profile is drawn in the XZ plane and revolved about Z. It starts at the
    disc's bottom center, walks out around the disc, follows the exponential
    down to 5 mm, and closes through the spike tip back along the axis.
    """
    disc_cone_profile = (
        cq.Workplane("XZ")
        .moveTo(0, 0)
        .lineTo(DISC_OUTER_RADIUS_MILLIMETERS, 0)
        .lineTo(DISC_OUTER_RADIUS_MILLIMETERS, DISC_THICKNESS_MILLIMETERS)
    )
    for segment_index in range(1, EXPONENTIAL_CONE_PROFILE_SEGMENT_COUNT + 1):
        height_fraction = segment_index / EXPONENTIAL_CONE_PROFILE_SEGMENT_COUNT
        profile_radius_millimeters = exponential_cone_radius_millimeters(height_fraction)
        profile_height_millimeters = (
            DISC_THICKNESS_MILLIMETERS + NARROWING_CONE_HEIGHT_MILLIMETERS * height_fraction
        )
        disc_cone_profile = disc_cone_profile.lineTo(
            profile_radius_millimeters,
            profile_height_millimeters,
        )
    revolved_disc_cone = (
        disc_cone_profile.lineTo(0, DISC_CONE_OVERALL_HEIGHT_MILLIMETERS).close().revolve(360)
    )
    return revolved_disc_cone.cut(make_arch_gap_join_cutouts(arch_count))


def uncut_disc_cone_volume_cubic_millimeters() -> float:
    """Analytic volume of the disc, exponential cone, and spike before cutouts.

    The cone term is the integral of pi * r(z)^2 for r(z) = R0 * exp(k * z / H).
    Cutouts remove a little of this, so the finished solid is slightly smaller.
    """
    disc_volume_cubic_millimeters = (
        math.pi
        * DISC_OUTER_RADIUS_MILLIMETERS
        * DISC_OUTER_RADIUS_MILLIMETERS
        * DISC_THICKNESS_MILLIMETERS
    )
    narrowing_cone_volume_cubic_millimeters = (
        math.pi
        * NARROWING_CONE_HEIGHT_MILLIMETERS
        * (
            NARROWING_CONE_END_RADIUS_MILLIMETERS * NARROWING_CONE_END_RADIUS_MILLIMETERS
            - DISC_OUTER_RADIUS_MILLIMETERS * DISC_OUTER_RADIUS_MILLIMETERS
        )
        / (2.0 * EXPONENTIAL_CONE_DECAY_CONSTANT)
    )
    spike_volume_cubic_millimeters = (
        math.pi
        * NARROWING_CONE_END_RADIUS_MILLIMETERS
        * NARROWING_CONE_END_RADIUS_MILLIMETERS
        * SPIKE_HEIGHT_MILLIMETERS
        / 3.0
    )
    return (
        disc_volume_cubic_millimeters
        + narrowing_cone_volume_cubic_millimeters
        + spike_volume_cubic_millimeters
    )


def verify_cylinder_seats_in_cutouts(
    disc_cone: cq.Workplane,
    hollow_cylinder: cq.Workplane,
    arch_count: int,
) -> None:
    """Place the cylinder on the cone and check that only the feet occupy cutouts.

    Same orientation for both parts: arch 0 is on +X. The cylinder is shifted
    up so its bottom face lies on the cutout floors. A real clash shows up as
    intersection volume. Each gap center must be empty cone and solid cylinder,
    and each arch center must be the reverse so the cone still passes through
    the arch opening.
    """
    disc_cone_shape = disc_cone.val()
    seated_hollow_cylinder = hollow_cylinder.translate(
        (0.0, 0.0, seated_cylinder_bottom_height_millimeters())
    )
    seated_overlap = disc_cone.intersect(seated_hollow_cylinder)
    seated_overlap_volume_cubic_millimeters = sum(
        overlap_solid.Volume() for overlap_solid in seated_overlap.solids().vals()
    )
    if seated_overlap_volume_cubic_millimeters > MAXIMUM_SEATED_OVERLAP_VOLUME_CUBIC_MILLIMETERS:
        raise SystemExit(
            "seated cylinder intersects the cone by "
            f"{seated_overlap_volume_cubic_millimeters:.2f} mm^3"
        )

    cutout_floor_height_millimeters = seated_cylinder_bottom_height_millimeters()
    seated_cylinder_shape = seated_hollow_cylinder.val()
    sample_height_inside_cutout_millimeters = (
        cutout_floor_height_millimeters + JOIN_CUTOUT_DEPTH_MILLIMETERS / 2.0
    )
    for arch_index in range(arch_count):
        arch_gap_angle_radians = math.radians(360.0 * (arch_index + 0.5) / arch_count)
        arch_station_angle_radians = math.radians(360.0 * arch_index / arch_count)
        arch_gap_foot_sample = (
            CYLINDER_MIDWALL_RADIUS_MILLIMETERS * math.cos(arch_gap_angle_radians),
            CYLINDER_MIDWALL_RADIUS_MILLIMETERS * math.sin(arch_gap_angle_radians),
            sample_height_inside_cutout_millimeters,
        )
        arch_opening_sample = (
            CYLINDER_MIDWALL_RADIUS_MILLIMETERS * math.cos(arch_station_angle_radians),
            CYLINDER_MIDWALL_RADIUS_MILLIMETERS * math.sin(arch_station_angle_radians),
            sample_height_inside_cutout_millimeters,
        )
        foot_is_inside_cone = point_is_inside_solid(disc_cone_shape, *arch_gap_foot_sample)
        foot_is_inside_cylinder = point_is_inside_solid(
            seated_cylinder_shape, *arch_gap_foot_sample
        )
        if foot_is_inside_cone or not foot_is_inside_cylinder:
            raise SystemExit(f"arch gap {arch_index} does not sit in its cutout")
        arch_opening_is_inside_cone = point_is_inside_solid(disc_cone_shape, *arch_opening_sample)
        arch_opening_is_inside_cylinder = point_is_inside_solid(
            seated_cylinder_shape, *arch_opening_sample
        )
        if not arch_opening_is_inside_cone or arch_opening_is_inside_cylinder:
            raise SystemExit(f"cone does not pass through arch {arch_index}")
        # Just under the pocket the cone must still be solid, so the foot has a floor.
        cutout_floor_sample = (
            CYLINDER_MIDWALL_RADIUS_MILLIMETERS * math.cos(arch_gap_angle_radians),
            CYLINDER_MIDWALL_RADIUS_MILLIMETERS * math.sin(arch_gap_angle_radians),
            cutout_floor_height_millimeters - 0.2,
        )
        if not point_is_inside_solid(disc_cone_shape, *cutout_floor_sample):
            raise SystemExit(f"cutout {arch_index} has no floor for the arch gap to land on")


def verify_disc_cone(disc_cone: cq.Workplane, arch_count: int) -> None:
    """Check the disc, the exponential taper, the spike, and the join cutouts."""
    disc_cone_shape = disc_cone.val()
    if not disc_cone_shape.isValid():
        raise SystemExit("disc cone solid is not valid")

    bounding_box = disc_cone_shape.BoundingBox()
    expected_axis_limits_millimeters = {
        "x": (bounding_box.xmin, bounding_box.xmax, -DISC_OUTER_RADIUS_MILLIMETERS, DISC_OUTER_RADIUS_MILLIMETERS),
        "y": (bounding_box.ymin, bounding_box.ymax, -DISC_OUTER_RADIUS_MILLIMETERS, DISC_OUTER_RADIUS_MILLIMETERS),
        "z": (bounding_box.zmin, bounding_box.zmax, 0.0, DISC_CONE_OVERALL_HEIGHT_MILLIMETERS),
    }
    for axis_name, (
        actual_low_millimeters,
        actual_high_millimeters,
        expected_low_millimeters,
        expected_high_millimeters,
    ) in expected_axis_limits_millimeters.items():
        low_error_millimeters = abs(actual_low_millimeters - expected_low_millimeters)
        high_error_millimeters = abs(actual_high_millimeters - expected_high_millimeters)
        if (
            low_error_millimeters > BOUNDING_BOX_TOLERANCE_MILLIMETERS
            or high_error_millimeters > BOUNDING_BOX_TOLERANCE_MILLIMETERS
        ):
            raise SystemExit(
                f"disc cone bounding box {axis_name} is "
                f"{actual_low_millimeters:.3f}..{actual_high_millimeters:.3f}"
            )

    narrowing_cone_mid_height_millimeters = (
        DISC_THICKNESS_MILLIMETERS + NARROWING_CONE_HEIGHT_MILLIMETERS / 2.0
    )
    narrowing_cone_mid_radius_millimeters = exponential_cone_radius_millimeters(0.5)
    # A straight taper would still contain this point. The exponential does not.
    straight_taper_mid_radius_millimeters = (
        DISC_OUTER_RADIUS_MILLIMETERS + NARROWING_CONE_END_RADIUS_MILLIMETERS
    ) / 2.0
    spike_mid_height_millimeters = (
        DISC_THICKNESS_MILLIMETERS
        + NARROWING_CONE_HEIGHT_MILLIMETERS
        + SPIKE_HEIGHT_MILLIMETERS / 2.0
    )
    spike_mid_radius_millimeters = NARROWING_CONE_END_RADIUS_MILLIMETERS / 2.0
    inside_outside_samples = [
        (0.0, 0.0, DISC_THICKNESS_MILLIMETERS / 2.0, True, "disc center"),
        (DISC_OUTER_RADIUS_MILLIMETERS - 0.5, 0.0, DISC_THICKNESS_MILLIMETERS / 2.0, True, "disc rim"),
        (DISC_OUTER_RADIUS_MILLIMETERS - 0.5, 0.0, DISC_THICKNESS_MILLIMETERS + 1.0, False, "above the disc rim"),
        (
            narrowing_cone_mid_radius_millimeters - 1.0,
            0.0,
            narrowing_cone_mid_height_millimeters,
            True,
            "inside the cone",
        ),
        (
            narrowing_cone_mid_radius_millimeters + 1.0,
            0.0,
            narrowing_cone_mid_height_millimeters,
            False,
            "outside the cone",
        ),
        (
            straight_taper_mid_radius_millimeters,
            0.0,
            narrowing_cone_mid_height_millimeters,
            False,
            "outside the exponential profile",
        ),
        (
            NARROWING_CONE_END_RADIUS_MILLIMETERS - 0.3,
            0.0,
            DISC_THICKNESS_MILLIMETERS + NARROWING_CONE_HEIGHT_MILLIMETERS - 0.3,
            True,
            "cone top",
        ),
        (spike_mid_radius_millimeters - 0.2, 0.0, spike_mid_height_millimeters, True, "inside the spike"),
        (spike_mid_radius_millimeters + 0.4, 0.0, spike_mid_height_millimeters, False, "outside the spike"),
        (0.0, 0.0, DISC_CONE_OVERALL_HEIGHT_MILLIMETERS + 1.0, False, "above the tip"),
    ]
    cutout_sample_height_millimeters = (
        seated_cylinder_bottom_height_millimeters() + JOIN_CUTOUT_DEPTH_MILLIMETERS / 2.0
    )
    first_arch_gap_angle_radians = math.pi / arch_count
    arch_gap_x_unit = math.cos(first_arch_gap_angle_radians)
    arch_gap_y_unit = math.sin(first_arch_gap_angle_radians)
    inside_outside_samples.extend(
        [
            (
                CYLINDER_MIDWALL_RADIUS_MILLIMETERS * arch_gap_x_unit,
                CYLINDER_MIDWALL_RADIUS_MILLIMETERS * arch_gap_y_unit,
                cutout_sample_height_millimeters,
                False,
                "join cutout at the cylinder bore",
            ),
            (
                (CYLINDER_INNER_RADIUS_MILLIMETERS - JOIN_FIT_CLEARANCE_MILLIMETERS - 0.25)
                * arch_gap_x_unit,
                (CYLINDER_INNER_RADIUS_MILLIMETERS - JOIN_FIT_CLEARANCE_MILLIMETERS - 0.25)
                * arch_gap_y_unit,
                cutout_sample_height_millimeters,
                True,
                "inside the cutout's bore side",
            ),
            (
                CYLINDER_MIDWALL_RADIUS_MILLIMETERS * arch_gap_x_unit,
                CYLINDER_MIDWALL_RADIUS_MILLIMETERS * arch_gap_y_unit,
                seated_cylinder_bottom_height_millimeters() - 0.2,
                True,
                "cutout floor",
            ),
            (
                CYLINDER_MIDWALL_RADIUS_MILLIMETERS,
                0.0,
                cutout_sample_height_millimeters,
                True,
                "cone remains in the arch opening",
            ),
            (
                CYLINDER_MIDWALL_RADIUS_MILLIMETERS * arch_gap_x_unit,
                CYLINDER_MIDWALL_RADIUS_MILLIMETERS * arch_gap_y_unit,
                DISC_THICKNESS_MILLIMETERS + 0.75,
                True,
                "cone base is no longer notched",
            ),
        ]
    )
    for height_fraction in (0.25, 0.5, 0.75):
        profile_radius_millimeters = exponential_cone_radius_millimeters(height_fraction)
        profile_height_millimeters = (
            DISC_THICKNESS_MILLIMETERS + NARROWING_CONE_HEIGHT_MILLIMETERS * height_fraction
        )
        inside_outside_samples.append(
            (
                profile_radius_millimeters - 0.4,
                0.0,
                profile_height_millimeters,
                True,
                f"inside cone at {height_fraction:.2f}",
            )
        )
        inside_outside_samples.append(
            (
                profile_radius_millimeters + 0.4,
                0.0,
                profile_height_millimeters,
                False,
                f"outside cone at {height_fraction:.2f}",
            )
        )
    for (
        sample_x_millimeters,
        sample_y_millimeters,
        sample_z_millimeters,
        sample_should_be_inside,
        sample_label,
    ) in inside_outside_samples:
        sample_is_inside = point_is_inside_solid(
            disc_cone_shape,
            sample_x_millimeters,
            sample_y_millimeters,
            sample_z_millimeters,
        )
        if sample_is_inside != sample_should_be_inside:
            raise SystemExit(
                f"{sample_label}: expected inside={sample_should_be_inside}, got {sample_is_inside}"
            )

    uncut_volume_cubic_millimeters = uncut_disc_cone_volume_cubic_millimeters()
    finished_volume_cubic_millimeters = disc_cone_shape.Volume()
    # Cutouts remove material, but only a small fraction of the disc-cone.
    volume_is_plausible = (
        uncut_volume_cubic_millimeters * 0.9
        < finished_volume_cubic_millimeters
        < uncut_volume_cubic_millimeters
    )
    if not volume_is_plausible:
        raise SystemExit(
            f"disc cone volume {finished_volume_cubic_millimeters:.1f} mm^3, "
            f"uncut {uncut_volume_cubic_millimeters:.1f}"
        )


def export_shape_to_stl(shape: cq.Workplane, stl_path: Path) -> None:
    """Write an STL using the shared linear and angular deflections."""
    cq.exporters.export(
        shape,
        str(stl_path),
        tolerance=STL_LINEAR_DEFLECTION_MILLIMETERS,
        angularTolerance=STL_ANGULAR_DEFLECTION_RADIANS,
    )


def main() -> None:
    arch_count = choose_arch_count()
    hollow_cylinder, arch_count = build_hollow_cylinder()
    verify_hollow_cylinder(hollow_cylinder, arch_count)
    export_shape_to_stl(hollow_cylinder, HOLLOW_CYLINDER_STL_PATH)
    midwall_arch_gap_millimeters = clear_arch_gap_arc_length_millimeters(
        CYLINDER_MIDWALL_RADIUS_MILLIMETERS,
        arch_count,
    )
    print(f"wrote {HOLLOW_CYLINDER_STL_FILENAME}")
    print(f"outer diameter {CYLINDER_OUTER_DIAMETER_MILLIMETERS:.2f} mm")
    print(f"inner diameter {CYLINDER_INNER_DIAMETER_MILLIMETERS:.2f} mm")
    print(f"height {CYLINDER_HEIGHT_MILLIMETERS:.2f} mm")
    print(f"wall {CYLINDER_WALL_THICKNESS_MILLIMETERS:.2f} mm")
    print(f"arches {arch_count}, radius {ARCH_CUTOUT_RADIUS_MILLIMETERS:.2f} mm")
    print(
        f"gap at mid-wall {midwall_arch_gap_millimeters:.3f} mm "
        f"(requested {REQUESTED_ARCH_GAP_MILLIMETERS:.2f} mm)"
    )
    print(
        "gap at inner surface "
        f"{clear_arch_gap_arc_length_millimeters(CYLINDER_INNER_RADIUS_MILLIMETERS, arch_count):.3f} mm"
    )
    print(
        "gap at outer surface "
        f"{clear_arch_gap_arc_length_millimeters(CYLINDER_OUTER_RADIUS_MILLIMETERS, arch_count):.3f} mm"
    )
    print(f"volume {hollow_cylinder.val().Volume():.1f} mm^3")

    disc_cone = make_disc_cone(arch_count)
    verify_disc_cone(disc_cone, arch_count)
    verify_cylinder_seats_in_cutouts(disc_cone, hollow_cylinder, arch_count)
    export_shape_to_stl(disc_cone, DISC_CONE_STL_PATH)
    print(f"wrote {DISC_CONE_STL_FILENAME}")
    print(
        f"disc diameter {DISC_OUTER_DIAMETER_MILLIMETERS:.2f} mm, "
        f"thickness {DISC_THICKNESS_MILLIMETERS:.2f} mm"
    )
    print(
        f"exponential cone height {NARROWING_CONE_HEIGHT_MILLIMETERS:.2f} mm, "
        f"end diameter {NARROWING_CONE_END_DIAMETER_MILLIMETERS:.2f} mm"
    )
    print(
        "diameter at mid-cone "
        f"{2.0 * exponential_cone_radius_millimeters(0.5):.2f} mm"
    )
    print(
        f"spike height {SPIKE_HEIGHT_MILLIMETERS:.2f} mm, "
        f"base diameter {SPIKE_BASE_DIAMETER_MILLIMETERS:.2f} mm"
    )
    cutout_floor_height_above_disc_millimeters = (
        seated_cylinder_bottom_height_millimeters() - DISC_THICKNESS_MILLIMETERS
    )
    print(
        f"join cutouts {arch_count} for bore {CYLINDER_INNER_DIAMETER_MILLIMETERS:.2f} mm, "
        f"{cutout_floor_height_above_disc_millimeters:.2f} mm above the disc"
    )
    print(
        f"cutout depth {JOIN_CUTOUT_DEPTH_MILLIMETERS:.2f} mm, "
        f"clearance {JOIN_FIT_CLEARANCE_MILLIMETERS:.2f} mm per side"
    )
    print(f"overall height {DISC_CONE_OVERALL_HEIGHT_MILLIMETERS:.2f} mm")
    print(f"volume {disc_cone.val().Volume():.1f} mm^3")


if __name__ == "__main__":
    main()