from pathlib import Path
import math

import cadquery as cq


AXIS = cq.Vector(0.0, 1.0, 0.0)


def _point_on_xz(radius: float, angle_degrees: float, y: float) -> cq.Vector:
    angle = math.radians(angle_degrees)
    return cq.Vector(radius * math.cos(angle), y, radius * math.sin(angle))


def _cylinder(radius: float, y0: float, length: float, x: float = 0.0, z: float = 0.0):
    return cq.Solid.makeCylinder(radius, length, cq.Vector(x, y0, z), AXIS)


def _ring(outer_radius: float, inner_radius: float, y0: float, length: float):
    return _cylinder(outer_radius, y0, length).cut(_cylinder(inner_radius, y0, length))


def _conical_ring(
    outer_start: float,
    outer_end: float,
    inner_start: float,
    inner_end: float,
    y0: float,
    length: float,
):
    outer = cq.Solid.makeCone(outer_start, outer_end, length, cq.Vector(0.0, y0, 0.0), AXIS)
    inner = cq.Solid.makeCone(inner_start, inner_end, length, cq.Vector(0.0, y0, 0.0), AXIS)
    return outer.cut(inner)


def _sector_cut(angle_center: float, angle_width: float, y0: float, length: float, radius: float = 100.0):
    start = math.radians(angle_center - angle_width / 2.0)
    end = math.radians(angle_center + angle_width / 2.0)
    points = [
        (0.0, 0.0),
        (radius * math.cos(start), radius * math.sin(start)),
        (radius * math.cos(end), radius * math.sin(end)),
    ]
    plane = cq.Plane(origin=(0.0, y0, 0.0), xDir=(1.0, 0.0, 0.0), normal=(0.0, 1.0, 0.0))
    return cq.Workplane(plane).polyline(points).close().extrude(length).val()


def _annular_sector_cut(
    angle_center: float,
    angle_width: float,
    inner_radius: float,
    outer_radius: float,
    y0: float,
    length: float,
    profile_fillet_radius: float = 0.0,
    axial_fillet_radius: float = 0.0,
):
    start = math.radians(angle_center - angle_width / 2.0)
    middle = math.radians(angle_center)
    end = math.radians(angle_center + angle_width / 2.0)
    inner_start = (inner_radius * math.cos(start), inner_radius * math.sin(start))
    inner_middle = (inner_radius * math.cos(middle), inner_radius * math.sin(middle))
    inner_end = (inner_radius * math.cos(end), inner_radius * math.sin(end))
    outer_start = (outer_radius * math.cos(start), outer_radius * math.sin(start))
    outer_middle = (outer_radius * math.cos(middle), outer_radius * math.sin(middle))
    outer_end = (outer_radius * math.cos(end), outer_radius * math.sin(end))
    plane = cq.Plane(origin=(0.0, y0, 0.0), xDir=(1.0, 0.0, 0.0), normal=(0.0, 1.0, 0.0))
    cutter = (
        cq.Workplane(plane)
        .moveTo(*inner_start)
        .lineTo(*outer_start)
        .threePointArc(outer_middle, outer_end)
        .lineTo(*inner_end)
        .threePointArc(inner_middle, inner_start)
        .close()
        .extrude(length)
    )
    if profile_fillet_radius > 0.0:
        # Round the corners of the annular profile (fillet axes parallel to Y).
        cutter = cutter.edges("|Y").fillet(profile_fillet_radius)
    if axial_fillet_radius > 0.0:
        # Round only the straight radial edges at the two axial ends.  This
        # produces the R10/R5/R2 window run-outs present in the source STEP.
        end_edges = [
            edge
            for edge in cutter.edges("%LINE").vals()
            if abs(edge.Center().y - y0) < 1.0e-6
            or abs(edge.Center().y - (y0 + length)) < 1.0e-6
        ]
        cutter = cutter.newObject(end_edges).fillet(axial_fillet_radius)
    return cutter.val()


def _radial_bar(angle_degrees: float, y0: float, length: float, inner_radius: float, outer_radius: float, width: float):
    radial_length = outer_radius - inner_radius
    bar = (
        cq.Workplane("XY")
        .box(radial_length, length, width)
        .translate(((inner_radius + outer_radius) / 2.0, y0 + length / 2.0, 0.0))
    )
    return bar.rotate((0.0, 0.0, 0.0), (0.0, 1.0, 0.0), -angle_degrees).val()


def _rounded_radial_slot(angle_degrees: float, y0: float, length: float):
    """A 40 mm wide radial window with R2 corners in its axial section."""
    plane = cq.Plane(
        origin=(40.0, y0 + length / 2.0, 0.0),
        xDir=(0.0, 1.0, 0.0),
        normal=(1.0, 0.0, 0.0),
    )
    slot = (
        cq.Workplane(plane)
        .rect(length - 4.0, 36.0)
        .offset2D(2.0)
        .extrude(42.0)
        .rotate((0.0, 0.0, 0.0), (0.0, 1.0, 0.0), -angle_degrees)
    )
    return slot.val()


def _support_blend(angle_degrees: float, mirrored: bool):
    """One of the six R2 blends between the R64 and R29.2 support arcs."""
    sharp = (-62.80715975198, 12.298808230444)
    inner_tangent = (-62.1855700706093, 15.131254904777)
    lobe_tangent = (-63.6707499640005, 16.8512614004542)
    inner_center = (0.0, 0.0)
    blend_center = (-61.6824120316291, 17.0669284278203)
    lobe_center = (-34.6410161514, 20.0)

    def transform(point):
        x, z = point
        if mirrored:
            axis = math.radians(150.0)
            ux, uz = math.cos(axis), math.sin(axis)
            projection = x * ux + z * uz
            x, z = 2.0 * projection * ux - x, 2.0 * projection * uz - z
        angle = math.radians(angle_degrees)
        return (
            x * math.cos(angle) - z * math.sin(angle),
            x * math.sin(angle) + z * math.cos(angle),
        )

    def arc_middle(center, start, end, radius):
        ux = (start[0] - center[0]) / radius + (end[0] - center[0]) / radius
        uz = (start[1] - center[1]) / radius + (end[1] - center[1]) / radius
        scale = radius / math.hypot(ux, uz)
        return center[0] + ux * scale, center[1] + uz * scale

    points = [transform(point) for point in (sharp, inner_tangent, lobe_tangent)]
    centers = [transform(point) for point in (inner_center, blend_center, lobe_center)]
    plane = cq.Plane(origin=(0.0, 90.2, 0.0), xDir=(1.0, 0.0, 0.0), normal=(0.0, -1.0, 0.0))
    return (
        cq.Workplane(plane)
        .moveTo(*points[0])
        .threePointArc(arc_middle(centers[0], points[0], points[1], 64.0), points[1])
        .threePointArc(arc_middle(centers[1], points[1], points[2], 2.0), points[2])
        .threePointArc(arc_middle(centers[2], points[2], points[0], 29.2), points[0])
        .close()
        .extrude(10.0)
        .val()
    )


def build():
    # Front cylindrical shell and its smaller continuation.
    part = cq.Workplane(obj=_ring(79.4, 74.4, -10.0, 22.0))
    part = part.union(_ring(78.4, 74.4, 12.0, 43.2))

    # Nine radial R2 holes through the front lip, spaced every 40 degrees.
    for angle in (90.0, 50.0, 10.0, -30.0, -70.0, -110.0, -150.0, -190.0, -230.0):
        start = _point_on_xz(72.0, angle, -5.0)
        direction = _point_on_xz(1.0, angle, 0.0)
        direction = cq.Vector(direction.x, 0.0, direction.z)
        part = part.cut(cq.Solid.makeCylinder(2.0, 10.0, start, direction))

    # Six long ventilation windows, 40-degree openings with 20-degree ribs.
    for center in (-110.0, -50.0, 10.0, 70.0, 130.0, 190.0):
        part = part.cut(
            _annular_sector_cut(center, 40.0, 70.0, 82.0, 12.0, 40.0, axial_fillet_radius=10.0)
        )

    # Taper from the large front shell to the narrower middle shell.
    part = part.union(_conical_ring(78.4, 73.2, 74.4, 69.2, 55.2, 20.0))
    part = part.union(_ring(73.2, 60.0, 75.2, 5.0))

    # Three small axial holes in the first shoulder ring.
    for angle in (30.0, 150.0, 270.0):
        p = _point_on_xz(64.0, angle, 75.2)
        part = part.cut(cq.Solid.makeCylinder(2.0, 5.0, p, AXIS))

    # Three-lobed support plate: a disk with three large circular openings.
    support = cq.Workplane(obj=_ring(73.2, 64.0, 80.2, 10.0))
    for x, z in ((34.6410161514, 20.0), (-34.6410161514, 20.0), (0.0, -40.0)):
        support = support.cut(_cylinder(29.2, 80.2, 10.0, x, z))
    part = part.union(support)

    # Thin middle shell and its second, angularly shifted six-window pattern.
    part = part.union(_ring(73.2, 69.2, 90.2, 30.0))
    for center in (-150.0, -90.0, -30.0, 30.0, 90.0, 150.0):
        part = part.cut(
            _annular_sector_cut(center, 40.0, 65.0, 76.0, 90.2, 20.0, axial_fillet_radius=5.0)
        )

    # Rear shell is one continuous 2:7 taper.  Its inner cavity changes from
    # a steep R69.2 -> R63.6 cone to the matching shallow taper at Y=125.8.
    taper_outer = cq.Solid.makeCone(
        73.2, 64.91428571428571, 29.0, cq.Vector(0.0, 120.2, 0.0), AXIS
    )
    taper_inner_a = cq.Solid.makeCone(69.2, 63.6, 5.6, cq.Vector(0.0, 120.2, 0.0), AXIS)
    taper_inner_b = cq.Solid.makeCone(63.6, 59.2, 15.4, cq.Vector(0.0, 125.8, 0.0), AXIS)
    part = part.union(taper_outer.cut(taper_inner_a).cut(taper_inner_b))
    hub = cq.Workplane(obj=_ring(38.0, 26.0, 129.2, 4.0)).union(
        _ring(38.0, 30.0, 133.2, 16.0)
    )
    part = part.union(hub)

    # Stop 0.0001 mm before the coincident rear-plate face so OpenCascade
    # keeps the six ribs in one valid solid.
    for center in (30.0, 90.0, 150.0, 210.0, 270.0, 330.0):
        part = part.cut(_rounded_radial_slot(center, 131.2, 9.9999))

    # Central bore through the solid rear plate, then its six fan windows.
    part = part.cut(_cylinder(30.0, 141.2, 8.0))
    for center in (45.0, 105.0, 165.0, 225.0, 285.0, 345.0):
        cutter = _annular_sector_cut(
            center, 30.0, 38.0, 59.2, 141.2, 8.0, profile_fillet_radius=2.0
        )
        part = part.cut(cutter)

    # These holes pass through both the end plate and the previously added hub.
    for angle in (-90.0, -30.0, 30.0, 90.0, 150.0, 210.0):
        p = _point_on_xz(34.0, angle, 141.2)
        part = part.cut(cq.Solid.makeCylinder(1.81198, 8.0, p, AXIS))

    # Add the small support blends last so they cannot interfere with fusing
    # the face-connected rear shell and end plate.
    for mirrored in (False, True):
        for angle in (0.0, 120.0, 240.0):
            part = part.union(_support_blend(angle, mirrored))

    # Remove only the six 0.0001 mm run-out residuals after the complete part
    # is connected.  This reaches the exact Y=141.2 boundary without causing
    # the earlier coincident-face split into six solids.
    for center in (30.0, 90.0, 150.0, 210.0, 270.0, 330.0):
        full_slot = _rounded_radial_slot(center, 131.2, 10.0)
        stable_slot = _rounded_radial_slot(center, 131.2, 9.9999)
        part = part.cut(full_slot.cut(stable_slot))

    return part


if __name__ == "__main__":
    output = Path(__file__).resolve().with_name("外壳框架.step")
    cq.exporters.export(build(), str(output))
    print(output)