AR-Shipdesign/arshipdesign/parametric/wizard_cruiser.py

"""
Generador paramétrico — Casco de Desplazamiento / Crucero (carena redonda).

Genera una tabla de offsets para embarcaciones de desplazamiento:
  - Carena redondeada (round bilge)
  - Sección maestra llena (Cm 0.82–0.90)
  - Proa fina, popa de espejo o crucero
  - LCB ajustable

Parámetros típicos:
  Cb: 0.45 – 0.65
  Velocidad/Froude: Fn 0.20 – 0.35

Autor: Álvaro Romero  |  Sprint 2A — AR-ShipDesign
"""
from __future__ import annotations

import math

import numpy as np

from arshipdesign.core.hull import Hull
from arshipdesign.core.offsets import OffsetsTable


# ---------------------------------------------------------------------------
# Forma de sección — carena redonda tipo desplazamiento
# ---------------------------------------------------------------------------

def _round_bilge_section(z: float, T: float, y_wl: float, Cm: float) -> float:
    """
    Sección de carena redonda parametrizada por Cm.

    Usa la forma  y = y_wl · sin(π/2 · (z/T)ⁿ)  donde n se obtiene de Cm.

    - n < 1  →  sección llena (Cm alto), convexa hacia afuera
    - n = 1  →  arco de círculo (Cm ≈ 0.637)
    - n > 1  →  sección en V (Cm bajo), convexa hacia el eje

    Propiedades garantizadas:
        y(z=0) = 0    (quilla puntual)
        y(z=T) = y_wl (manga en flotación)
        dy/dz ≥ 0     (monótona, sin inflexión indeseada)
    """
    if y_wl < 1e-9 or T < 1e-9:
        return 0.0
    t = min(1.0, max(0.0, z / T))
    if t < 1e-12:
        return 0.0
    # Tabla calibrada Cm → n (verificada contra integración numérica)
    _CM = [0.45, 0.55, 0.65, 0.75, 0.82, 0.88, 0.94]
    _N  = [2.20, 1.60, 1.15, 0.80, 0.58, 0.42, 0.28]
    n = float(np.interp(Cm, _CM, _N))
    return float(y_wl * math.sin(math.pi / 2.0 * t ** n))


def make_displacement_hull(
    name: str = "Crucero de Desplazamiento",
    lpp: float = 12.0,
    beam: float = 3.80,
    draft: float = 1.40,
    depth: float = 2.20,
    cb: float = 0.55,
    lcb_frac: float = 0.52,        # fracción de Lpp desde AP
    cm: float = 0.86,              # coeficiente de cuaderna maestra
    n_stations: int = 21,
    n_waterlines: int = 11,
) -> Hull:
    """Genera un casco de desplazamiento de carena redonda.

    Parámetros
    ----------
    cb : float
        Coeficiente de bloque objetivo (0.45–0.65).
    lcb_frac : float
        Posición del LCB como fracción de Lpp desde AP (0.50–0.55).
    cm : float
        Coeficiente de cuaderna maestra (0.82–0.92).
    """
    x_sta = np.linspace(0.0, lpp, n_stations)
    z_wl  = np.linspace(0.0, draft, n_waterlines)
    xi    = (x_sta / lpp - 0.5) * 2.0   # ∈ [−1, 1], 0=midship

    # ── Plan form (semi-manga en flotación) ────────────────────────────
    # LCB desplazado del midship
    lcb_shift = 2.0 * (lcb_frac - 0.5)   # ∈ [−1, 1]
    f_plan = _displacement_plan_form(xi, cb, lcb_shift)

    data = np.zeros((n_stations, n_waterlines))

    for i in range(n_stations):
        y_wl = (beam / 2.0) * f_plan[i]

        # Cm varía a lo largo de la eslora: lleno en midship, más en V en extremos.
        # f_plan[i]=1 → midship → Cm=cm; f_plan[i]→0 → extremos → Cm≈0.52
        local_cm = float(np.clip(
            cm * (0.42 + 0.58 * f_plan[i] ** 0.40),
            0.52, cm
        ))

        for j, z in enumerate(z_wl):
            data[i, j] = _round_bilge_section(z, draft, y_wl, local_cm)

    data = np.clip(data, 0.0, None)

    offsets = OffsetsTable(
        x_stations=x_sta,
        z_waterlines=z_wl,
        data=data,
        station_labels=[f"S{i}" for i in range(n_stations)],
        lpp=lpp, beam=beam, draft=draft,
    )
    return Hull(
        name=name, lpp=lpp, beam=beam, depth=depth, draft=draft, offsets=offsets
    )


def _displacement_plan_form(
    xi: np.ndarray, cb: float, lcb_shift: float
) -> np.ndarray:
    """
    Plan form normalizada para carena de desplazamiento.

    Usa una distribución tipo seno modificado con el apex desplazado
    según la posición del LCB.

    xi ∈ [−1, 1], −1=AP, +1=FP.
    Retorna f ∈ [0, 1].
    """
    # Exponent de la plan form: más alto = más llena = Cb mayor
    n = _plan_exponent(cb)

    xi_shifted = xi - lcb_shift
    xi_shifted = np.clip(xi_shifted, -1.0, 1.0)

    # Plan form base: bell curve asimétrica
    # Parte de entrada (FP): más suave
    # Parte de salida (AP): ligeramente más llena
    f = np.zeros_like(xi)
    for k, x in enumerate(xi_shifted):
        if x <= 0:
            # Run (desde midship hacia AP): ligeramente más llena
            t = -x
            f[k] = max(0.0, 1.0 - (t ** (n * 0.90)))
        else:
            # Entry (desde midship hacia FP): estándar
            t = x
            f[k] = max(0.0, 1.0 - (t ** n))

    return f


def _plan_exponent(cb: float) -> float:
    """Mapa Cb → exponente de plan form."""
    cb_vals = [0.35, 0.45, 0.55, 0.65, 0.75]
    n_vals  = [1.2,  1.6,  2.0,  2.7,  3.5]
    return float(np.interp(cb, cb_vals, n_vals))