"""
Formed limb strip — after press brake operations on the flat blank.

Three forming stages:
  1. Belly curve: gentle reflex along the full working limb
  2. Siyah recurve: 70 deg forward bend over last ~120mm at each tip
  3. String retention: conical taper with rolled edges.
     The tip narrows to a cone in plan view.
     The belly-side edges roll upward and inward toward the centerline,
     forming retention lips that cradle the string loop.
     String drops in from above, self-centers on the cone taper
     under tension.

Center section (inside riser, +/-140mm) stays flat.

Run: blender --background --python ambi-recurve-limb-formed.py
"""

import bpy
import bmesh
import math


def clear_scene():
    bpy.ops.object.select_all(action='SELECT')
    bpy.ops.object.delete()


def make_material(name, color, metallic=0.0, roughness=0.5):
    mat = bpy.data.materials.new(name)
    mat.use_nodes = True
    bsdf = mat.node_tree.nodes["Principled BSDF"]
    bsdf.inputs["Base Color"].default_value = (*color, 1.0)
    bsdf.inputs["Metallic"].default_value = metallic
    bsdf.inputs["Roughness"].default_value = roughness
    return mat


# ---------------------------------------------------------------
# Taper — Gaussian bell curve
# ---------------------------------------------------------------

def limb_taper(d_mm):
    """(half_width_mm, thickness_mm) at distance d_mm from center."""
    d = abs(d_mm)
    w_peak, w_floor, w_sigma = 22.0, 10.0, 450.0
    t_peak, t_floor, t_sigma = 3.5, 1.4, 420.0
    hw = w_floor + (w_peak - w_floor) * math.exp(-d*d / (2*w_sigma*w_sigma))
    th = t_floor + (t_peak - t_floor) * math.exp(-d*d / (2*t_sigma*t_sigma))
    return (hw, th)


# ---------------------------------------------------------------
# Path (belly + siyah)
# ---------------------------------------------------------------

BELLY_START  = 140.0
BELLY_END    = 360.0
BELLY_RADIUS = 800.0
SIYAH_START  = BELLY_END
SIYAH_END    = 480.0
SIYAH_ANGLE  = math.radians(70)
SIYAH_ARC    = SIYAH_END - SIYAH_START
SIYAH_RADIUS = SIYAH_ARC / SIYAH_ANGLE


def limb_path(d_mm):
    """(y_mm, z_mm, angle_rad) along the neutral axis."""
    if d_mm <= BELLY_START:
        return (d_mm, 0.0, 0.0)

    belly_arc = BELLY_END - BELLY_START
    belly_theta_end = belly_arc / BELLY_RADIUS
    y_belly_end = BELLY_START + BELLY_RADIUS * math.sin(belly_theta_end)
    z_belly_end = BELLY_RADIUS * (1 - math.cos(belly_theta_end))

    if d_mm <= BELLY_END:
        a = (d_mm - BELLY_START) / BELLY_RADIUS
        return (BELLY_START + BELLY_RADIUS * math.sin(a),
                BELLY_RADIUS * (1 - math.cos(a)), a)

    base_a = belly_theta_end
    if d_mm <= SIYAH_END:
        a = (d_mm - SIYAH_START) / SIYAH_RADIUS
        total_a = base_a + a
        dy = SIYAH_RADIUS * (math.sin(total_a) - math.sin(base_a))
        dz = SIYAH_RADIUS * (math.cos(base_a) - math.cos(total_a))
        return (y_belly_end + dy, z_belly_end + dz, total_a)

    return limb_path(SIYAH_END)


def tip_state():
    return limb_path(SIYAH_END)


# ---------------------------------------------------------------
# Build main strip body (one arm)
# ---------------------------------------------------------------

def build_strip_body(bm, sign, steps=100):
    """Build one arm. Returns list of (v0,v1,v2,v3) rings."""
    dd = SIYAH_END / steps
    rings = []
    for i in range(steps + 1):
        d = i * dd
        hw_mm, th_mm = limb_taper(d)
        y_mm, z_mm, angle = limb_path(d)

        y = sign * y_mm / 1000.0
        z = z_mm / 1000.0
        hw = hw_mm / 1000.0
        th = th_mm / 1000.0

        ny = -math.sin(angle) * sign
        nz = math.cos(angle)

        v0 = bm.verts.new((-hw, y,           z))
        v1 = bm.verts.new(( hw, y,           z))
        v2 = bm.verts.new(( hw, y + ny * th, z + nz * th))
        v3 = bm.verts.new((-hw, y + ny * th, z + nz * th))
        rings.append((v0, v1, v2, v3))

    for i in range(len(rings) - 1):
        for j in range(4):
            jn = (j + 1) % 4
            bm.faces.new([rings[i][j], rings[i][jn],
                          rings[i+1][jn], rings[i+1][j]])
    return rings


# ---------------------------------------------------------------
# Build cone tip with rolled edges for one arm
# ---------------------------------------------------------------

def build_cone_tip(bm, sign):
    """
    Conical taper with S-curve edge profile for string retention.

    Edge profile (looking end-on at one edge):

        head (small, tight)
         __
        /  \    <-- string sits HERE, under the head
       |    |
        \   |
         \  |   <-- belly (big, sweeping arc) — clearance zone
          \ |       wire only touches the head, not the limb
           \|
            |   <-- strip body

    The S has:
      - Big belly: large radius arc sweeping forward, creating
        clearance between wire and limb surface
      - Small head: tight reverse curve hooking back over the belly,
        forming the spring lip that holds the string
      - The string wire touches ONLY the head
      - The head acts as a leaf spring: flexes under load

    The area between belly and limb body is the narrow gap that
    keeps the wire clear. Cone taper self-centers the loop.
    """
    ty_mm, tz_mm, t_angle = tip_state()
    hw_mm, th_mm = limb_taper(SIYAH_END)

    cos_a = math.cos(t_angle)
    sin_a = math.sin(t_angle)
    tg_y = cos_a * sign
    tg_z = sin_a
    nm_y = -sin_a * sign
    nm_z = cos_a
    fw_y = sin_a * sign
    fw_z = -cos_a

    by = sign * ty_mm / 1000.0
    bz = tz_mm / 1000.0

    def world_pt(lat_mm, along_mm, up_mm, fwd_mm=0):
        x = lat_mm / 1000.0
        y = (by + (along_mm / 1000.0) * tg_y
                + (up_mm / 1000.0) * nm_y
                + (fwd_mm / 1000.0) * fw_y)
        z = (bz + (along_mm / 1000.0) * tg_z
                + (up_mm / 1000.0) * nm_z
                + (fwd_mm / 1000.0) * fw_z)
        return bm.verts.new((x, y, z))

    cone_length = 25.0
    steps_along = 24
    tip_hw = 2.0

    # S-curve parameters
    belly_r = 3.0     # big belly radius (mm) — the large sweeping arc
    head_r = 1.0      # small head radius (mm) — the tight hook-back

    # S-curve profile for one edge (in the fwd/up plane):
    #
    # Start at strip corner: (fwd=0, up=0)
    # Belly arc: large radius, sweeps forward and up
    #   center at (0, belly_r), sweeps from -90 deg to ~+90 deg
    #   -> endpoint at (belly_r, belly_r) roughly
    #   but we only go partway: about 160 deg of arc
    # Head arc: small radius, reverses back
    #   tangent-continuous with belly, hooks back over
    #   -> creates the lip overhang

    # Points per side of the S-curve
    belly_pts = 8    # points along the big belly arc
    head_pts = 5     # points along the small head arc
    pts_per_side = belly_pts + head_pts
    # Total ring: left S + bottom center + right S + 3 back spine
    pts_per_ring = pts_per_side * 2 + 1 + 3

    all_rings = []

    for i in range(steps_along + 1):
        frac = i / steps_along
        d = frac * cone_length
        frac_s = frac * frac * (3 - 2 * frac)

        local_hw = hw_mm + (tip_hw - hw_mm) * frac_s

        # How developed the S-curve is: 0 at base (flat), 1 at tip (full S)
        s_dev = frac_s

        # Belly arc: sweep angle increases from 0 to 160 deg
        belly_sweep = math.radians(160) * s_dev
        # Head arc: sweep angle increases from 0 to 140 deg
        head_sweep = math.radians(140) * s_dev

        ring_verts = []

        def s_curve_points(side):
            """
            Generate S-curve profile points for one edge.
            side: -1 for left, +1 for right.
            Returns list of (fwd_mm, up_mm) tuples.
            """
            pts = []

            # --- Belly arc ---
            # Center of belly arc: at the edge corner, offset inward by belly_r
            # Arc starts at (fwd=0, up=0) and sweeps forward and up
            # Starting angle: -90 deg (pointing down = at the corner)
            # Sweep: belly_sweep degrees counterclockwise
            belly_cx = 0.0  # fwd center
            belly_cz = belly_r  # up center (above the corner)

            for k in range(belly_pts):
                t = k / (belly_pts - 1)  # 0 to 1
                a = -math.pi / 2 + belly_sweep * t
                fwd = belly_cx + belly_r * math.cos(a)
                up = belly_cz + belly_r * math.sin(a)
                pts.append((fwd, up))

            if belly_sweep > 0.01 and head_sweep > 0.01:
                # --- Head arc ---
                # Tangent-continuous with belly end.
                # At the belly endpoint, the tangent direction is
                # perpendicular to the radius at that point.
                belly_end_angle = -math.pi / 2 + belly_sweep
                belly_end_fwd = belly_cx + belly_r * math.cos(belly_end_angle)
                belly_end_up = belly_cz + belly_r * math.sin(belly_end_angle)

                # The head arc center is offset from the belly endpoint
                # in the opposite direction of the belly radius
                # (to create the S reversal)
                # Belly radius direction at endpoint:
                br_fwd = math.cos(belly_end_angle)
                br_up = math.sin(belly_end_angle)
                # Head center: step INWARD from belly end by head_r
                # in the radius direction
                head_cx = belly_end_fwd + head_r * br_fwd
                head_cz = belly_end_up + head_r * br_up

                # Head arc starts at the belly endpoint
                # and sweeps in the OPPOSITE rotational direction
                # Starting angle for head: belly_end_angle + pi
                # (pointing back toward the belly endpoint)
                head_start = belly_end_angle + math.pi

                for k in range(head_pts):
                    t = k / (head_pts - 1)
                    a = head_start - head_sweep * t
                    fwd = head_cx + head_r * math.cos(a)
                    up = head_cz + head_r * math.sin(a)
                    pts.append((fwd, up))

            return pts

        def s_curve_points_with_lip(side, local_hw_val):
            """
            Generate S-curve + lip. The lip is a lateral flare on the
            INSIDE face of the head — the metal is stretched sideways
            (toward strip center) to create a smooth bearing shelf
            where the wire sits. Prevents the edge from cutting
            into the wire rope.
            """
            base_pts = s_curve_points(side)
            result = []

            n_belly = belly_pts
            for idx, (fwd, up) in enumerate(base_pts):
                lat = side * local_hw_val  # default lateral position

                # For the head section (after belly points),
                # the inner face gets a lateral lip/shelf
                if idx >= n_belly:
                    head_idx = idx - n_belly
                    head_t = head_idx / max(head_pts - 1, 1)
                    # Lip grows from 0 to ~0.8mm inward on the
                    # inner face of the head. Eased.
                    lip_t = head_t * head_t * (3 - 2 * head_t)
                    lip_inward = 0.8 * lip_t  # mm toward center
                    lat = side * (local_hw_val - lip_inward)

                result.append((lat, fwd, up))

            return result

            return pts

        # --- Left S-curve with lip ---
        left_pts = s_curve_points_with_lip(-1, local_hw)
        for lat, fwd, up in left_pts:
            ring_verts.append(world_pt(lat, d, up, fwd))

        # --- Bottom center ---
        ring_verts.append(world_pt(0, d, 0, 0))

        # --- Right S-curve with lip (reversed for continuous ring) ---
        right_pts = s_curve_points_with_lip(+1, local_hw)
        for lat, fwd, up in reversed(right_pts):
            ring_verts.append(world_pt(lat, d, up, fwd))

        # --- Back spine ---
        ring_verts.append(world_pt(local_hw, d, th_mm, 0))
        ring_verts.append(world_pt(0, d, th_mm, 0))
        ring_verts.append(world_pt(-local_hw, d, th_mm, 0))

        all_rings.append(ring_verts)

    # Faces between consecutive rings
    for i in range(len(all_rings) - 1):
        n = len(all_rings[i])
        if len(all_rings[i+1]) != n:
            continue
        for j in range(n):
            jn = (j + 1) % n
            bm.faces.new([all_rings[i][j], all_rings[i][jn],
                          all_rings[i+1][jn], all_rings[i+1][j]])

    # Cap the tip
    tip_ring = all_rings[-1]
    center_v = world_pt(0, cone_length, th_mm / 2, 0)
    n = len(tip_ring)
    for j in range(n):
        jn = (j + 1) % n
        bm.faces.new([tip_ring[j], tip_ring[jn], center_v])


# ---------------------------------------------------------------
# Main
# ---------------------------------------------------------------

def main():
    clear_scene()
    steel = make_material("316L_QPQ", (0.08, 0.08, 0.08), metallic=0.8, roughness=0.4)

    bm = bmesh.new()

    upper = build_strip_body(bm, +1)
    lower = build_strip_body(bm, -1)

    # bridge center
    u, l = upper[0], lower[0]
    for j in range(4):
        jn = (j + 1) % 4
        bm.faces.new([u[j], u[jn], l[jn], l[j]])

    # cone tips with rolled edges
    build_cone_tip(bm, +1)
    build_cone_tip(bm, -1)

    mesh = bpy.data.meshes.new("Limb_Formed")
    bm.to_mesh(mesh)
    bm.free()

    obj = bpy.data.objects.new("Limb_Strip_Formed", mesh)
    bpy.context.scene.collection.objects.link(obj)
    obj.data.materials.append(steel)

    for f in obj.data.polygons:
        f.use_smooth = True

    mod = obj.modifiers.new("Subsurf", 'SUBSURF')
    mod.levels = 2
    mod.render_levels = 3

    # bolt holes
    for offset_mm in [-30, 0, 30]:
        bpy.ops.mesh.primitive_cylinder_add(
            radius=0.00325, depth=0.008,
            location=(0, offset_mm / 1000.0, 0.00175))
        hole = bpy.context.active_object
        hole.name = f"BoltHole_{offset_mm}"
        bmod = obj.modifiers.new(f"Hole_{offset_mm}", 'BOOLEAN')
        bmod.operation = 'DIFFERENCE'
        bmod.object = hole
        bpy.context.view_layer.objects.active = obj
        bpy.ops.object.modifier_apply(modifier=f"Hole_{offset_mm}")
        bpy.data.objects.remove(hole)

    bpy.ops.object.select_all(action='DESELECT')
    obj.select_set(True)
    bpy.context.view_layer.objects.active = obj

    print("Done: formed limb strip with S-curve tip retention.")
    print("  Taper: Gaussian bell curve")
    print("  1. Belly curve: R800mm, 140-360mm")
    print("  2. Siyah recurve: 70 deg, 360-480mm")
    print("  3. Cone tip with S-curve edges: 480-505mm")
    print("     - Plan view: narrows 20mm -> 4mm (cone)")
    print("     - Edge profile: S-curve (big belly + small head)")
    print("     - Belly (R3mm): sweeps forward, clearance zone")
    print("     - Head (R1mm): hooks back, spring lip holds string")
    print("     - Wire touches ONLY the head, clear of limb body")


if __name__ == "__main__":
    main()