"""
gz_integrator.py — Motor de cálculo de la curva GZ de estabilidad estática.

Dos métodos:
  1. compute_gz_wall_sided  — Fórmula de pared lateral (rápida, exacta hasta ~30°).
  2. compute_gz_direct      — Integración directa 3-D con recorte de polígono
                              por Sutherland-Hodgman (precisa a grandes ángulos).

Convenciones
------------
  * Todos los ángulos de escora φ en grados en la API pública;
    en radianes internamente.
  * Sistema de coordenadas del casco (upright):
      y = manga (+ estribor), z = altura (0 = quilla, + hacia arriba).
  * Sistema de coordenadas del mundo (heeled):
      Eje z_w vertical apuntando arriba; eje y_w horizontal apuntando
      a estribor mundo.
  * La escora φ es positiva a estribor.

Fórmula GZ
----------
    GZ = y_B_world − KG · sin(φ)
  donde y_B_world es la coordenada horizontal (eje y_w) del centro de
  carena B en el sistema del mundo.

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

import math
from dataclasses import dataclass, field
from typing import Optional

import numpy as np
from scipy.integrate import simpson as _scipy_simpson
from scipy.optimize import brentq


# ---------------------------------------------------------------------------
# Dataclasses públicos
# ---------------------------------------------------------------------------

@dataclass
class GZPoint:
    """Un punto de la curva GZ."""
    phi_deg: float
    gz: float


@dataclass
class GZCurve:
    """Curva GZ completa con estadísticas derivadas.

    Atributos
    ---------
    hull_name : str
    draft : float  [m]
    kg : float     [m] — Altura del centro de gravedad sobre la quilla.
    gm : float     [m] — Altura metacéntrica transversal GM = KMT - KG.
    bmt : float    [m] — Radio metacéntrico transversal.
    points : list[GZPoint]
    area_0_30 : float   [m·rad]
    area_30_40 : float  [m·rad]
    area_0_40 : float   [m·rad]
    gz_max : float      [m]
    phi_gz_max : float  [°]
    avs : float         [°] — Ángulo de estabilidad nula (Angle of Vanishing Stability).
    """
    hull_name: str
    draft: float
    kg: float
    gm: float
    bmt: float
    points: list[GZPoint] = field(default_factory=list)
    area_0_30: float = 0.0
    area_30_40: float = 0.0
    area_0_40: float = 0.0
    gz_max: float = 0.0
    phi_gz_max: float = 0.0
    avs: float = 90.0

    # ------------------------------------------------------------------
    # Propiedades de acceso rápido a arrays
    # ------------------------------------------------------------------

    @property
    def angles_deg(self) -> np.ndarray:
        return np.array([p.phi_deg for p in self.points], dtype=float)

    @property
    def gz_values(self) -> np.ndarray:
        return np.array([p.gz for p in self.points], dtype=float)

    # ------------------------------------------------------------------
    # Cálculo de áreas, gz_max y AVS
    # ------------------------------------------------------------------

    def _compute_areas(self) -> None:
        """Integra las áreas bajo la curva GZ y calcula gz_max, phi_gz_max, avs.

        Usa scipy.integrate.simpson sobre los puntos disponibles.
        Las áreas se calculan en [m·rad] (ángulos convertidos a radianes).
        """
        phi = self.angles_deg
        gz  = self.gz_values

        if len(phi) < 2:
            return

        # Convertir ángulos a radianes para la integración
        phi_rad = np.deg2rad(phi)

        # ── Área 0–30° ──────────────────────────────────────────────
        mask_030 = phi <= 30.0
        if mask_030.sum() >= 2:
            self.area_0_30 = float(abs(_scipy_simpson(
                gz[mask_030], x=phi_rad[mask_030]
            )))
        else:
            self.area_0_30 = 0.0

        # ── Área 30–40° ─────────────────────────────────────────────
        # Incluir el punto exacto en 30° (interpolado si no existe)
        gz30 = float(np.interp(30.0, phi, gz))
        gz40 = float(np.interp(40.0, phi, gz))

        mask_3040 = (phi >= 30.0) & (phi <= 40.0)
        pts_3040_phi = np.concatenate([[30.0], phi[mask_3040], [40.0]])
        pts_3040_gz  = np.concatenate([[gz30], gz[mask_3040], [gz40]])

        # Eliminar duplicados
        _, unique_idx = np.unique(pts_3040_phi, return_index=True)
        pts_3040_phi = pts_3040_phi[unique_idx]
        pts_3040_gz  = pts_3040_gz[unique_idx]

        if len(pts_3040_phi) >= 2:
            self.area_30_40 = float(abs(_scipy_simpson(
                pts_3040_gz, x=np.deg2rad(pts_3040_phi)
            )))
        else:
            self.area_30_40 = 0.0

        # ── Área 0–40° ──────────────────────────────────────────────
        mask_040 = phi <= 40.0
        pts_040_phi = np.concatenate([phi[mask_040], [40.0]])
        pts_040_gz  = np.concatenate([gz[mask_040], [gz40]])
        _, unique_idx2 = np.unique(pts_040_phi, return_index=True)
        pts_040_phi = pts_040_phi[unique_idx2]
        pts_040_gz  = pts_040_gz[unique_idx2]

        if len(pts_040_phi) >= 2:
            self.area_0_40 = float(abs(_scipy_simpson(
                pts_040_gz, x=np.deg2rad(pts_040_phi)
            )))
        else:
            self.area_0_40 = 0.0

        # ── GZ máximo y su ángulo ────────────────────────────────────
        idx_max = int(np.argmax(gz))
        self.gz_max = float(gz[idx_max])
        self.phi_gz_max = float(phi[idx_max])

        # ── AVS: primer cruce de GZ = 0 después del pico ────────────
        # Buscar a partir del índice del máximo
        avs = 90.0
        for i in range(idx_max, len(gz) - 1):
            if gz[i] > 0.0 and gz[i + 1] <= 0.0:
                # Interpolación lineal para encontrar el cruce exacto
                phi1, gz1 = phi[i],     gz[i]
                phi2, gz2 = phi[i + 1], gz[i + 1]
                if abs(gz2 - gz1) > 1e-12:
                    avs = phi1 + (phi2 - phi1) * (-gz1) / (gz2 - gz1)
                else:
                    avs = float(phi1)
                break
        self.avs = avs


# ---------------------------------------------------------------------------
# Método 1: Fórmula de pared lateral (Wall-Sided)
# ---------------------------------------------------------------------------

def compute_gz_wall_sided(
    hull,
    draft: float,
    kg: Optional[float] = None,
    angles_deg: Optional[np.ndarray] = None,
    rho: float = 1025.0,
) -> GZCurve:
    """Calcula la curva GZ por la fórmula de pared lateral.

    GZ = sin(φ) · (GM + 0.5 · BM · tan²(φ))

    Esta fórmula es exacta para cascos de paredes verticales (wall-sided)
    a cualquier ángulo; para cascos con flare proporciona buena aproximación
    hasta ~30° y subestima a ángulos mayores.

    Parameters
    ----------
    hull : Hull
    draft : float  [m]
    kg : float, optional  — Si None, usa hull.depth * 0.55
    angles_deg : array_like, optional — Default: 0..90 en pasos de 1°
    rho : float — Densidad del agua [kg/m³]

    Returns
    -------
    GZCurve
    """
    from arshipdesign.hydrostatics.upright import compute_upright

    hydro = compute_upright(hull, draft, rho=rho)

    kg_val = float(hull.depth * 0.55 if kg is None else kg)
    gm     = hydro.kmt - kg_val
    bmt    = hydro.bmt

    if angles_deg is None:
        phi_arr = np.linspace(0.0, 90.0, 91)
    else:
        phi_arr = np.asarray(angles_deg, dtype=float)

    points: list[GZPoint] = []
    for phi_deg in phi_arr:
        phi_rad = math.radians(phi_deg)
        sin_phi = math.sin(phi_rad)
        tan_phi = math.tan(phi_rad) if abs(phi_rad) < math.radians(89.9) else math.tan(math.radians(89.9))
        gz = sin_phi * (gm + 0.5 * bmt * tan_phi ** 2)
        points.append(GZPoint(phi_deg=float(phi_deg), gz=float(gz)))

    curve = GZCurve(
        hull_name=hull.name,
        draft=float(draft),
        kg=kg_val,
        gm=gm,
        bmt=bmt,
        points=points,
    )
    curve._compute_areas()
    return curve


# ---------------------------------------------------------------------------
# Método 2: Integración directa (Direct)
# ---------------------------------------------------------------------------

def _clip_polygon_below_z(
    poly: list[tuple[float, float]],
    z_wl: float,
) -> list[tuple[float, float]]:
    """Sutherland-Hodgman clip: mantiene región z ≤ z_wl.

    Parameters
    ----------
    poly : list of (y, z)
    z_wl : float — línea de corte (cota de flotación en el sistema girado)

    Returns
    -------
    list of (y, z) — polígono recortado (puede ser vacío)
    """
    if not poly:
        return []

    output = list(poly)

    def _inside(pt: tuple[float, float]) -> bool:
        return pt[1] <= z_wl

    def _intersect(
        a: tuple[float, float], b: tuple[float, float]
    ) -> tuple[float, float]:
        """Intersección de segmento (a→b) con z = z_wl."""
        ya, za = a
        yb, zb = b
        dz = zb - za
        if abs(dz) < 1e-14:
            return ((ya + yb) * 0.5, z_wl)
        t = (z_wl - za) / dz
        return (ya + t * (yb - ya), z_wl)

    # Plano de recorte: z ≤ z_wl (mantener interior)
    input_list = output
    output = []

    n = len(input_list)
    for i in range(n):
        curr = input_list[i]
        prev = input_list[i - 1]
        if _inside(curr):
            if not _inside(prev):
                output.append(_intersect(prev, curr))
            output.append(curr)
        elif _inside(prev):
            output.append(_intersect(prev, curr))

    return output


def _polygon_area_centroid(
    poly: list[tuple[float, float]],
) -> tuple[float, float, float]:
    """Área y centroide de un polígono cerrado por la fórmula de Shoelace.

    Returns
    -------
    area : float  — Área absoluta (siempre ≥ 0).
    yc : float    — Coordenada y del centroide.
    zc : float    — Coordenada z del centroide.
    """
    n = len(poly)
    if n < 3:
        if n == 0:
            return 0.0, 0.0, 0.0
        y0, z0 = poly[0]
        return 0.0, float(y0), float(z0)

    # Shoelace
    ys = np.array([p[0] for p in poly], dtype=float)
    zs = np.array([p[1] for p in poly], dtype=float)

    # Área con signo
    cross = ys[:-1] * zs[1:] - ys[1:] * zs[:-1]
    # Cerrar el polígono
    cross_close = ys[-1] * zs[0] - ys[0] * zs[-1]
    total_cross = np.sum(cross) + cross_close
    area_signed = 0.5 * total_cross
    area = abs(area_signed)

    if area < 1e-15:
        return 0.0, float(np.mean(ys)), float(np.mean(zs))

    # Centroides
    cx_terms = (ys[:-1] + ys[1:]) * cross
    cz_terms = (zs[:-1] + zs[1:]) * cross
    cx_close = (ys[-1] + ys[0]) * cross_close
    cz_close = (zs[-1] + zs[0]) * cross_close

    yc = (np.sum(cx_terms) + cx_close) / (6.0 * area_signed)
    zc = (np.sum(cz_terms) + cz_close) / (6.0 * area_signed)

    return area, float(yc), float(zc)


def _build_section_polygon_extended(
    half_breadths: np.ndarray,
    z_positions: np.ndarray,
    z_extend: float,
) -> list[tuple[float, float]]:
    """Construye el polígono extendido de una sección transversal.

    Para capturar correctamente el cálculo de GZ a grandes ángulos,
    el polígono se extiende con paredes verticales (wall-sided) desde
    la línea de flotación upright hasta z_extend.  Esto garantiza que
    la integral de Sutherland-Hodgman captura el cuño de inmersión /
    emersión al escorar.

    Geometría:
      - Estribor: de quilla (j=0) hacia arriba hasta el tope de la sección,
        luego pared vertical hasta z_extend.
      - Babor: pared vertical de z_extend hacia abajo, luego espejo de estribor
        hasta la quilla.
      - Resultado: 2n+2 vértices (n puntos de sección × 2 bandas + 2 extensiones).

    Parameters
    ----------
    half_breadths : ndarray, shape (n,)
    z_positions : ndarray, shape (n,)  — creciente desde la quilla
    z_extend : float                   — altura máxima del polígono (m)

    Returns
    -------
    list of (y, z) — polígono cerrado en sentido antihorario
    """
    n = len(z_positions)
    y_top = float(half_breadths[-1])
    z_top = float(z_positions[-1])  # normalmente = draft

    poly: list[tuple[float, float]] = []
    # Estribor: de quilla (j=0) hacia arriba
    for j in range(n):
        poly.append((float(half_breadths[j]), float(z_positions[j])))
    # Extensión estribor hasta z_extend (wall-sided)
    if z_extend > z_top + 1e-9:
        poly.append((y_top, z_extend))
    # Extensión babor desde z_extend hacia abajo (wall-sided, mirror)
    if z_extend > z_top + 1e-9:
        poly.append((-y_top, z_extend))
    # Babor: de tope hacia abajo (mirror)
    for j in reversed(range(n)):
        poly.append((-float(half_breadths[j]), float(z_positions[j])))
    return poly


def _compute_polygon_volume(
    upright_polys: list[list[tuple[float, float]]],
    dx: np.ndarray,
    z_wl: float,
) -> float:
    """Suma del volumen de los polígonos recortados en z ≤ z_wl."""
    vol = 0.0
    for i, poly in enumerate(upright_polys):
        clipped = _clip_polygon_below_z(poly, z_wl)
        if clipped:
            a, _, _ = _polygon_area_centroid(clipped)
            vol += a * dx[i]
    return vol


def compute_gz_direct(
    hull,
    draft: float,
    kg: Optional[float] = None,
    angles_deg: Optional[np.ndarray] = None,
    rho: float = 1025.0,
) -> GZCurve:
    """Calcula la curva GZ por integración directa de secciones.

    Para cada ángulo de escora φ (escora a estribor, φ > 0):
      1. Construye el polígono extendido (wall-sided por encima del calado)
         de cada sección para capturar los cuños de inmersión/emersión.
      2. Aplica la rotación de casco-a-mundo (giro en sentido horario=escora
         a estribor):
             y_w =  y·cos(φ) + z·sin(φ)
             z_w = -y·sin(φ) + z·cos(φ)
      3. Determina z_wl (posición de la flotación en el mundo) por bisección
         (brentq) tal que el volumen sumergido = V_target.
      4. Calcula y_B_world = centroide horizontal del volumen recortado.
      5. GZ = y_B_world − KG·sin(φ)

    La rotación utilizada es la convención de estabilidad naval estándar:
    +φ = escora a estribor → el lado estribor (y_body > 0) baja.

    Parameters
    ----------
    hull : Hull
    draft : float  [m]
    kg : float, optional  — Si None, usa hull.depth * 0.55
    angles_deg : array_like, optional — Default: 0..90 en pasos de 1°
    rho : float

    Returns
    -------
    GZCurve
    """
    from arshipdesign.hydrostatics.upright import compute_upright

    hydro  = compute_upright(hull, draft, rho=rho)
    bmt    = hydro.bmt
    kmt    = hydro.kmt
    kg_val = float(hull.depth * 0.55 if kg is None else kg)
    gm     = kmt - kg_val

    if angles_deg is None:
        phi_arr = np.linspace(0.0, 90.0, 91)
    else:
        phi_arr = np.asarray(angles_deg, dtype=float)

    # ── Preparar secciones y pesos de integración trapezoidal ────────
    sections = hull.offsets.to_sections()
    x_sta    = np.array([s.x for s in sections], dtype=float)
    n_sta    = len(x_sta)

    dx = np.zeros(n_sta, dtype=float)
    if n_sta >= 2:
        dx[0]  = (x_sta[1]  - x_sta[0])  * 0.5
        dx[-1] = (x_sta[-1] - x_sta[-2]) * 0.5
        for i in range(1, n_sta - 1):
            dx[i] = (x_sta[i + 1] - x_sta[i - 1]) * 0.5

    # ── Polígonos extendidos en posición upright ──────────────────────
    # Se extienden con paredes verticales hasta z_extend para capturar
    # el cuño de inmersión a grandes ángulos de escora.
    # z_extend debe ser suficientemente mayor que draft + beam/2 * sin(90°).
    z_extend = float(draft) + float(hull.beam) * 1.5 + float(hull.depth) * 0.5
    upright_polys: list[list[tuple[float, float]]] = []
    for sec in sections:
        poly = _build_section_polygon_extended(
            sec.half_breadths, sec.z_positions, z_extend
        )
        upright_polys.append(poly)

    # ── V_target: volumen coherente con los polígonos ─────────────────
    # Usar el volumen que los polígonos dan al cortarlos en z ≤ draft,
    # para garantizar coherencia numérica entre phi=0 y phi>0.
    V_target = _compute_polygon_volume(upright_polys, dx, float(draft))

    # Límites conservadores para la bisección
    z_bisect_high = z_extend * 0.95  # justo por debajo del tope extendido

    points: list[GZPoint] = []

    for phi_deg in phi_arr:
        phi_rad = math.radians(float(phi_deg))
        cos_phi = math.cos(phi_rad)
        sin_phi = math.sin(phi_rad)

        if abs(phi_rad) < 1e-9:
            # φ = 0: GZ = 0 exacto por simetría
            points.append(GZPoint(phi_deg=float(phi_deg), gz=0.0))
            continue

        # ── Rotación de cuerpo a mundo para escora a estribor ────────
        # Convención naval: +φ → estribor baja (giro horario en plano y-z)
        #   y_w =  y·cos(φ) + z·sin(φ)
        #   z_w = -y·sin(φ) + z·cos(φ)
        rotated_polys: list[list[tuple[float, float]]] = []
        for poly in upright_polys:
            world_poly = [
                ( y * cos_phi + z * sin_phi,
                 -y * sin_phi + z * cos_phi)
                for y, z in poly
            ]
            rotated_polys.append(world_poly)

        # ── Bisección para z_wl en el sistema del mundo ───────────────
        def _vol_minus_target(z_wl: float) -> float:
            vol = 0.0
            for i in range(n_sta):
                clipped = _clip_polygon_below_z(rotated_polys[i], z_wl)
                if clipped:
                    a, _, _ = _polygon_area_centroid(clipped)
                    vol += a * dx[i]
            return vol - V_target

        z_b_low  = 1e-4
        z_b_high = z_bisect_high

        f_low  = _vol_minus_target(z_b_low)
        f_high = _vol_minus_target(z_b_high)

        if f_low * f_high > 0:
            # Fallback: fórmula wall-sided
            phi_clamp = min(phi_rad, math.radians(89.0))
            gz_ws = sin_phi * (gm + 0.5 * bmt * math.tan(phi_clamp) ** 2)
            points.append(GZPoint(phi_deg=float(phi_deg), gz=float(gz_ws)))
            continue

        try:
            z_wl = brentq(_vol_minus_target, z_b_low, z_b_high,
                          xtol=1e-5, maxiter=150)
        except ValueError:
            phi_clamp = min(phi_rad, math.radians(89.0))
            gz_ws = sin_phi * (gm + 0.5 * bmt * math.tan(phi_clamp) ** 2)
            points.append(GZPoint(phi_deg=float(phi_deg), gz=float(gz_ws)))
            continue

        # ── Centroide horizontal de carena en el mundo ────────────────
        vol_total = 0.0
        moment_y  = 0.0

        for i in range(n_sta):
            clipped = _clip_polygon_below_z(rotated_polys[i], z_wl)
            if clipped:
                a, yc, _ = _polygon_area_centroid(clipped)
                vol_total += a * dx[i]
                moment_y  += a * yc * dx[i]

        if vol_total > 1e-12:
            y_B_world = moment_y / vol_total
        else:
            y_B_world = 0.0

        # GZ = y_B_world - y_G_world
        # G = (y_body=0, z_body=KG); en mundo (CW): y_G_w = 0*cos + KG*sin = KG*sin
        gz = y_B_world - kg_val * sin_phi
        points.append(GZPoint(phi_deg=float(phi_deg), gz=float(gz)))

    curve = GZCurve(
        hull_name=hull.name,
        draft=float(draft),
        kg=kg_val,
        gm=gm,
        bmt=bmt,
        points=points,
    )
    curve._compute_areas()
    return curve