# TEST_11 改訂（球体筐体実験）用: 評価板の矩形輪郭と球を一体化した筐体メッシュ生成．
#
# 球単体では矩形の評価板と断面円の形が合わず，四隅の貫通か縁の隙間のどちらかが生じる．
# Fusion 360 のロフトに相当する「矩形 → 円」の遷移スカートをスクリプトで生成し，
# 球のお椀部分と 1 つの watertight ソリッドに縫合することで，評価板がフタとして
# 隙間なく嵌まる筐体を作る．
"""球＋ロフト一体型筐体（評価板輪郭に整合する閉キャビティ）STL を生成する．

構成（上から下へ，単位 mm）:
    1. リム矩形: 評価板輪郭から rim_inset だけ内側の矩形．z=rim_z（板の実体に
       食い込む高さ）に置き，板とのオーバーラップで隙間なくシールする
    2. ロフト部: 矩形 → 球断面円（z=cut_z）への極座標 Hermite ブレンド．
       球側の端点勾配を球面の接線に一致させ滑らかに接続する
    3. 球部: 半径 radius・中心 (0,0,center_z) の球の z ≤ cut_z 部分（お椀）

内面と外面（thickness だけオフセット）を生成し，リム帯で縫合した閉ソリッドを出力する．

使用例:
    uv run python scripts/generate_sphere_loft_housing.py \
        --output configs/meshes/球ロフト_r100.stl
"""

from __future__ import annotations

import argparse
from pathlib import Path

import numpy as np
import trimesh


def _rect_radius(theta: np.ndarray, half_x: float, half_y: float) -> np.ndarray:
    """角度 theta ごとの矩形境界までの半径（極座標表現）を返す．"""
    cos_t = np.abs(np.cos(theta))
    sin_t = np.abs(np.sin(theta))
    with np.errstate(divide="ignore"):
        rx = np.where(cos_t > 1e-12, half_x / cos_t, np.inf)
        ry = np.where(sin_t > 1e-12, half_y / sin_t, np.inf)
    return np.minimum(rx, ry)


# blend="radial" のリム接線の張り出し強さ（弦長に対する比）．
# 大きいほど天井が板の外へ大きく張り出し，装置が扁平になる．
RADIAL_TANGENT_SCALE = 0.6


def _loft_ring(
    u: float,
    rect_r: np.ndarray,
    cut_r: float,
    cut_z: float,
    rim_z: float,
    cut_d: float,
    radius: float,
    blend: str,
) -> tuple[np.ndarray, np.ndarray]:
    """ロフト部のパラメータ u における 1 リングの (半径, Z) を返す．

    "cubic" / "quintic" は z を等間隔に進めながら半径 ρ(z) を決める形なので，
    リムでの傾きは dρ/dz = 0（＝板に垂直に立ち上がる）しか表現できない．これが
    板直下に「喉」を作り，半径が大きいほどくびれて見える原因になる．
    "radial" は (ρ, z) の 2 次元パラメトリック Hermite に切り替え，リムの接線を
    **半径方向（水平）**に取ることで喉を無くす（接続点側は従来どおり球の接線に一致）．

    Args:
        u: 0（リム矩形）→ 1（球との接続円）のパラメータ．
        rect_r: リム矩形の角度ごとの半径．
        cut_r: 接続円の半径．
        cut_z: 接続平面の Z．
        rim_z: リムの Z．
        cut_d: 接続平面の球中心からの高さ（cut_z - center_z）．
        radius: 球内半径．
        blend: "cubic" / "quintic" / "radial"．

    Returns:
        (半径, Z) の組．いずれも角度ごとの配列．
    """
    span = cut_z - rim_z
    if blend == "radial":
        # 弦長に比例した接線ベクトルを両端に置く 2 次元 3 次 Hermite
        chord = np.hypot(cut_r - rect_r, span)
        scale = RADIAL_TANGENT_SCALE * chord
        h00 = 2 * u**3 - 3 * u**2 + 1
        h10 = u**3 - 2 * u**2 + u
        h01 = -2 * u**3 + 3 * u**2
        h11 = u**3 - u**2
        rho = h00 * rect_r + h10 * scale + h01 * cut_r + h11 * scale * (cut_d / radius)
        z = h00 * rim_z + h01 * cut_z + h11 * scale * (-cut_r / radius)
        return rho, z
    m1 = (-cut_d / cut_r) * span
    m2 = (-(radius**2) / cut_r**3) * span**2
    return _blend_radius(u, rect_r, cut_r, m1, m2, blend), np.full_like(
        rect_r, rim_z + u * span
    )


def _blend_radius(
    u: float,
    rect_r: np.ndarray,
    cut_r: float,
    m1: float,
    m2: float,
    blend: str,
) -> np.ndarray:
    """ロフト部のパラメータ u における半径を Hermite ブレンドで返す．

    Args:
        u: 0（リム矩形）→ 1（球との接続円）のパラメータ．
        rect_r: リム矩形の角度ごとの半径．
        cut_r: 接続円の半径．
        m1: 接続円における dρ/du（球面の接線に一致させる値）．
        m2: 接続円における d²ρ/du²（球面の曲率に一致させる値．quintic のみ使用）．
        blend: "cubic"（傾きのみ一致）または "quintic"（傾きと曲率を一致）．

    Returns:
        角度ごとの半径配列．

    Raises:
        ValueError: 未対応の blend が指定された場合．
    """
    if blend == "cubic":
        return (2 * u**3 - 3 * u**2 + 1) * rect_r + (-2 * u**3 + 3 * u**2) * cut_r + (
            u**3 - u**2
        ) * m1
    if blend == "quintic":
        # 5 次 Hermite: p(0)=rect_r, p'(0)=p''(0)=0, p(1)=cut_r, p'(1)=m1, p''(1)=m2
        # 端点で曲率まで球に一致するので，接合部に曲率の段差（見た目の折れ）が出ない
        h0 = 1 - 10 * u**3 + 15 * u**4 - 6 * u**5
        h3 = 10 * u**3 - 15 * u**4 + 6 * u**5
        h4 = -4 * u**3 + 7 * u**4 - 3 * u**5
        h5 = 0.5 * (u**3 - 2 * u**4 + u**5)
        return h0 * rect_r + h3 * cut_r + h4 * m1 + h5 * m2
    msg = f"未対応の blend です: {blend!r}（cubic / quintic のいずれか）"
    raise ValueError(msg)


def _build_surface(
    radius: float,
    center_z: float,
    half_x: float,
    half_y: float,
    rim_z: float,
    cut_z: float,
    n_theta: int,
    n_loft: int,
    n_sphere: int,
    blend: str = "cubic",
) -> tuple[np.ndarray, np.ndarray]:
    """リム矩形 → ロフト → 球お椀の 1 枚曲面（頂点とリング構造）を生成する．

    Returns:
        (vertices, rings): vertices は (N, 3)，rings は (n_ring, n_theta) の頂点 index．
        リングは上（リム）から下（極直前）の順．末尾に極頂点を 1 つ追加している．
    """
    theta = np.linspace(0.0, 2.0 * np.pi, n_theta, endpoint=False)
    rect_r = _rect_radius(theta, half_x, half_y)

    cut_d = cut_z - center_z
    cut_r = float(np.sqrt(radius**2 - cut_d**2))
    rings_list: list[np.ndarray] = []
    z_list: list[np.ndarray] = []

    # ロフト部（リム u=0 を含み，カット面 u=1 は球部の先頭リングに任せる）
    for u in np.linspace(0.0, 1.0, n_loft, endpoint=False):
        rho, z_ring = _loft_ring(u, rect_r, cut_r, cut_z, rim_z, cut_d, radius, blend)
        rings_list.append(rho)
        z_list.append(z_ring)

    # 球部（カット面 → 極直前）．角度パラメータで均等分割
    phi_cut = float(np.arccos(cut_d / radius))
    for phi in np.linspace(phi_cut, np.pi, n_sphere, endpoint=False):
        rings_list.append(np.full(n_theta, radius * np.sin(phi)))
        z_list.append(np.full(n_theta, center_z + radius * np.cos(phi)))

    n_ring = len(rings_list)
    vertices = np.empty((n_ring * n_theta + 1, 3), dtype=np.float64)
    rings = np.arange(n_ring * n_theta).reshape(n_ring, n_theta)
    for i, (rho, z) in enumerate(zip(rings_list, z_list, strict=True)):
        vertices[rings[i], 0] = rho * np.cos(theta)
        vertices[rings[i], 1] = rho * np.sin(theta)
        vertices[rings[i], 2] = z
    vertices[-1] = (0.0, 0.0, center_z - radius)  # 極頂点
    return vertices, rings


def _stitch_rings(rings: np.ndarray, pole_index: int) -> np.ndarray:
    """隣接リング間を四角形 2 分割で三角形化し，最下リングは極頂点へファン接続する．"""
    n_ring, n_theta = rings.shape
    faces: list[tuple[int, int, int]] = []
    for i in range(n_ring - 1):
        a, b = rings[i], rings[i + 1]
        for j in range(n_theta):
            k = (j + 1) % n_theta
            faces.append((a[j], b[j], b[k]))
            faces.append((a[j], b[k], a[k]))
    bottom = rings[-1]
    for j in range(n_theta):
        k = (j + 1) % n_theta
        faces.append((bottom[j], pole_index, bottom[k]))
    return np.asarray(faces, dtype=np.int64)


def build_loft_housing(
    radius: float,
    thickness: float,
    center_z: float,
    board_half_x: float,
    board_half_y: float,
    rim_inset: float,
    rim_z: float,
    cut_z: float,
    n_theta: int = 256,
    n_loft: int = 48,
    n_sphere: int = 160,
    blend: str = "cubic",
) -> trimesh.Trimesh:
    """球＋ロフト一体型筐体の閉ソリッド（内面＋外面＋リム帯）を生成する．

    Args:
        radius: 球内半径（mm）．
        thickness: 殻厚（mm）．外面は内面を径方向におおむね thickness オフセットした形．
        center_z: 球中心の Z 座標（mm）．X, Y は 0 固定．
        board_half_x: 評価板輪郭の X 半幅（mm）．
        board_half_y: 評価板輪郭の Y 半幅（mm）．
        rim_inset: リム矩形を板輪郭から内側へ寄せる量（mm）．板とのオーバーラップ代．
        rim_z: リム矩形の Z 座標（mm）．板の実体に食い込む高さに置いてシールする．
        cut_z: ロフトと球の接続平面の Z 座標（mm）．
        n_theta: 周方向分割数．
        n_loft: ロフト部の軸方向リング数．
        n_sphere: 球部の軸方向リング数．
        blend: ロフトの形状．"cubic" は球の接線だけを一致させる（既定・後方互換）．
            "quintic" は曲率まで一致させ，接合部の折れを無くす．いずれもリムでは
            板に垂直に立ち上がるため板直下に喉ができる．"radial" はリムの接線を
            半径方向（水平）に取り，喉を無くす．

    Returns:
        watertight な筐体ソリッドメッシュ（mm 単位）．

    Raises:
        ValueError: 生成メッシュが watertight にならなかった場合．
    """
    hx_in = board_half_x - rim_inset
    hy_in = board_half_y - rim_inset
    v_in, rings_in = _build_surface(
        radius, center_z, hx_in, hy_in, rim_z, cut_z, n_theta, n_loft, n_sphere, blend
    )
    v_out, rings_out = _build_surface(
        radius + thickness,
        center_z,
        hx_in + thickness,
        hy_in + thickness,
        rim_z,
        cut_z,
        n_theta,
        n_loft,
        n_sphere,
        blend,
    )

    f_in = _stitch_rings(rings_in, len(v_in) - 1)
    f_out = _stitch_rings(rings_out, len(v_out) - 1)

    offset = len(v_in)
    vertices = np.vstack([v_in, v_out])
    # リム帯: 内面上端リングと外面上端リングを接続して閉じる
    top_in, top_out = rings_in[0], rings_out[0] + offset
    band: list[tuple[int, int, int]] = []
    for j in range(n_theta):
        k = (j + 1) % n_theta
        band.append((top_in[j], top_out[j], top_out[k]))
        band.append((top_in[j], top_out[k], top_in[k]))
    faces = np.vstack([f_in, f_out + offset, np.asarray(band, dtype=np.int64)])

    mesh = trimesh.Trimesh(vertices=vertices, faces=faces, process=False)
    mesh.fix_normals()  # 閉ソリッドとして法線を外向きに統一（内面は自動的にキャビティ向き）
    if not mesh.is_watertight:
        raise ValueError("生成メッシュが watertight になりませんでした（縫合パラメータを確認）")
    return mesh


def _parse_args(argv: list[str] | None = None) -> argparse.Namespace:
    """コマンドライン引数をパースして返す．"""
    parser = argparse.ArgumentParser(
        prog="generate_sphere_loft_housing",
        description="評価板輪郭と球を一体化した筐体 STL（mm 単位）を生成する．",
    )
    parser.add_argument("--radius", type=float, default=100.0, help="球内半径（mm）．既定 100")
    parser.add_argument("--thickness", type=float, default=3.0, help="殻厚（mm）．既定 3")
    parser.add_argument(
        "--center-z", type=float, default=-92.0, help="球中心の Z 座標（mm）．既定 -92"
    )
    parser.add_argument(
        "--board-half", type=float, nargs=2, default=(49.5, 51.0), metavar=("HX", "HY"),
        help="評価板輪郭の半幅 X Y（mm）．既定 49.5 51（99×102mm 板）",
    )
    parser.add_argument(
        "--rim-inset",
        type=float,
        default=1.5,
        help="リムを板輪郭から内側へ寄せる量（mm）．既定 1.5",
    )
    parser.add_argument(
        "--rim-z", type=float, default=0.0,
        help="リムの Z 座標（mm）．既定 0（板裏面 -4.5 より 4.5mm 食い込みシール）",
    )
    parser.add_argument(
        "--cut-z", type=float, default=-40.0, help="ロフト⇄球の接続平面 Z（mm）．既定 -40"
    )
    parser.add_argument(
        "--blend", choices=("cubic", "quintic", "radial"), default="cubic",
        help="ロフトの形状．quintic は接合部で曲率まで球に一致させる．"
        "radial はリムで半径方向（水平）に出して板直下の喉を無くす．既定 cubic",
    )
    parser.add_argument(
        "--output",
        type=Path,
        default=Path("configs/meshes/球ロフト_r100.stl"),
        help="出力 STL パス．既定 configs/meshes/球ロフト_r100.stl",
    )
    return parser.parse_args(argv)


def main(argv: list[str] | None = None) -> int:
    """筐体を生成して STL 出力し，検証情報を表示する．"""
    args = _parse_args(argv)
    mesh = build_loft_housing(
        radius=args.radius,
        thickness=args.thickness,
        center_z=args.center_z,
        board_half_x=args.board_half[0],
        board_half_y=args.board_half[1],
        rim_inset=args.rim_inset,
        rim_z=args.rim_z,
        cut_z=args.cut_z,
        blend=args.blend,
    )
    args.output.parent.mkdir(parents=True, exist_ok=True)
    mesh.export(args.output)
    print(f"出力: {args.output}")
    print(f"faces: {len(mesh.faces)}, watertight: {mesh.is_watertight}")
    print(f"bounds (mm): {mesh.bounds.round(2).tolist()}")
    print(f"体積 (mm^3): {mesh.volume:.0f}（正値なら法線外向きで正常）")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())