AR-Shipdesign/arshipdesign/hydrostatics/upright.py

"""
UprightHydrostatics — hidrostáticos en condición vertical (quilla recta).

Calcula el conjunto completo de variables hidrostáticas para un calado
dado: volumen, desplazamiento, plano de flotación, metacentros, TPC, MCT
y los cuatro coeficientes de forma.

Conforme con IACS Rec.34 §4.3 — verificación analítica y §6 — trazabilidad.

Autor: Álvaro Romero
Módulo 2 — AR-ShipDesign
"""
from __future__ import annotations

from dataclasses import dataclass

import numpy as np

from arshipdesign.hydrostatics.integrator import (
    integrate,
    waterplane_strips,
    section_areas_and_centroids,
)


# ---------------------------------------------------------------------------
# Dataclass de resultado
# ---------------------------------------------------------------------------

@dataclass
class UprightHydrostatics:
    """Conjunto completo de hidrostáticos upright a un calado dado.

    Todos los valores referidos a:
      - Origen longitudinal: AP (x = 0)
      - Origen vertical: quilla (z = 0)
      - Plano de crujía: eje de simetría

    Atributos
    ---------
    draft : float
        Calado T [m].
    volume : float
        Volumen de desplazamiento V [m³].
    displacement : float
        Desplazamiento Δ [t]  (agua salada ρ = 1025 kg/m³ por defecto).
    awp : float
        Área del plano de flotación Awp [m²].
    lcb : float
        Centro longitudinal de carena desde AP LCB [m].
    lcf : float
        Centro longitudinal de flotación desde AP LCF [m].
    kb : float
        Altura del centro de carena sobre la quilla KB [m].
    it : float
        Segundo momento transversal del plano de flotación IT [m⁴].
    il : float
        Segundo momento longitudinal del plano de flotación IL [m⁴].
    bmt : float
        Radio metacéntrico transversal BM_T = IT / V [m].
    bml : float
        Radio metacéntrico longitudinal BM_L = IL / V [m].
    kmt : float
        Altura del metacentro transversal KM_T = KB + BM_T [m].
    kml : float
        Altura del metacentro longitudinal KM_L = KB + BM_L [m].
    tpc : float
        Toneladas por centímetro de inmersión TPC [t/cm].
    mct : float
        Momento para cambiar asiento 1 cm MCT [t·m/cm].
    cb : float
        Coeficiente de bloque Cb  [-].
    cw : float
        Coeficiente de plano de flotación Cw [-].
    cm : float
        Coeficiente de cuaderna maestra Cm [-].
    cp : float
        Coeficiente prismático Cp [-].
    """

    draft:        float
    volume:       float
    displacement: float
    awp:          float
    lcb:          float
    lcf:          float
    kb:           float
    it:           float
    il:           float
    bmt:          float
    bml:          float
    kmt:          float
    kml:          float
    tpc:          float
    mct:          float
    cb:           float
    cw:           float
    cm:           float
    cp:           float


# ---------------------------------------------------------------------------
# Función de cálculo
# ---------------------------------------------------------------------------

def compute_upright(
    hull,
    draft: float,
    rho:   float        = 1025.0,
    kg:    float | None = None,
) -> UprightHydrostatics:
    """Calcula todos los hidrostáticos upright para *hull* al calado *draft*.

    El cálculo se realiza en una sola pasada sobre las secciones, reutilizando
    los arrays intermedios para evitar redundancia.

    Parameters
    ----------
    hull : Hull
        Casco de referencia (objeto ``arshipdesign.core.hull.Hull``).
    draft : float
        Calado de cálculo T [m].  Si T ≤ 0, retorna ceros.
    rho : float
        Densidad del agua [kg/m³].  Default 1025 (agua salada).
    kg : float | None
        Altura del centro de gravedad KG [m].  Si None se estima
        como ``hull.depth × 0.55`` (buque en rosca, conservador).

    Returns
    -------
    UprightHydrostatics
        Todos los hidrostáticos al calado *draft*.
    """
    T = float(draft)

    # Caso degenerado (calado nulo o negativo)
    if T <= 1e-6:
        return _zero_hydrostatics(T)

    # ----------------------------------------------------------------
    # 1. Secciones transversales → volumen, LCB, KB
    # ----------------------------------------------------------------
    sections = hull.offsets.to_sections()
    x_s, areas, cz = section_areas_and_centroids(sections, T)

    vol   = abs(integrate(areas,         x_s))
    delta = vol * rho / 1000.0

    if vol > 1e-12:
        lcb = integrate(areas * x_s, x_s) / vol
        kb  = integrate(areas * cz,  x_s) / vol
    else:
        lcb = hull.lpp / 2.0
        kb  = T / 2.0

    # ----------------------------------------------------------------
    # 2. Plano de flotación → Awp, LCF, IT, IL
    # ----------------------------------------------------------------
    x_wl, y_wl = waterplane_strips(hull.offsets, T)
    strip = 2.0 * y_wl                          # ancho total a cada x

    awp = abs(integrate(strip, x_wl))

    if awp > 1e-12:
        lcf = integrate(strip * x_wl, x_wl) / awp
    else:
        lcf = hull.lpp / 2.0

    # IT = (2/3) · ∫ y³ dx   (Rawson & Tupper §3.2)
    it = abs(integrate((2.0 / 3.0) * y_wl ** 3, x_wl))

    # IL = ∫ 2y · (x − LCF)² dx
    il = abs(integrate(strip * (x_wl - lcf) ** 2, x_wl))

    # ----------------------------------------------------------------
    # 3. Radios metacéntricos
    # ----------------------------------------------------------------
    bmt = it / vol if vol > 1e-12 else 0.0
    bml = il / vol if vol > 1e-12 else 0.0
    kmt = kb + bmt
    kml = kb + bml

    # ----------------------------------------------------------------
    # 4. TPC y MCT
    # ----------------------------------------------------------------
    tpc = awp * rho / 100_000.0

    kg_val = hull.depth * 0.55 if kg is None else float(kg)
    gml    = max(kb + bml - kg_val, 0.0)
    mct    = delta * gml / (100.0 * hull.lpp) if hull.lpp > 1e-12 else 0.0

    # ----------------------------------------------------------------
    # 5. Coeficientes de forma
    # ----------------------------------------------------------------
    cb = vol / (hull.lpp * hull.beam * T) if (hull.lpp * hull.beam * T) > 1e-12 else 0.0
    cw = awp / (hull.lpp * hull.beam)     if (hull.lpp * hull.beam)      > 1e-12 else 0.0

    # Área de cuaderna maestra: interpolar en x_mid
    x_mid = hull.lpp / 2.0
    am    = float(np.interp(x_mid, x_s, areas))
    cm    = am / (hull.beam * T)      if (hull.beam * T)      > 1e-12 else 0.0
    cp    = vol / (am * hull.lpp)     if (am * hull.lpp)      > 1e-12 else 0.0

    return UprightHydrostatics(
        draft=T, volume=vol, displacement=delta,
        awp=awp, lcb=lcb, lcf=lcf, kb=kb,
        it=it,   il=il,
        bmt=bmt, bml=bml, kmt=kmt, kml=kml,
        tpc=tpc, mct=mct,
        cb=cb, cw=cw, cm=cm, cp=cp,
    )


# ---------------------------------------------------------------------------
# Auxiliar privado
# ---------------------------------------------------------------------------

def _zero_hydrostatics(draft: float) -> UprightHydrostatics:
    """Devuelve un UprightHydrostatics con todos los valores en cero."""
    return UprightHydrostatics(
        draft=draft, volume=0.0, displacement=0.0,
        awp=0.0, lcb=0.0, lcf=0.0, kb=0.0,
        it=0.0, il=0.0,
        bmt=0.0, bml=0.0, kmt=0.0, kml=0.0,
        tpc=0.0, mct=0.0,
        cb=0.0, cw=0.0, cm=0.0, cp=0.0,
    )