alro65
0f85935fc8
- gz_integrator.py: GZCurve, GZPoint, compute_gz_wall_sided (fórmula pared lateral), compute_gz_direct (integración Sutherland-Hodgman) - imo_is2008.py: IMOCriterion, IMOResult, check_imo_is2008 — 6 criterios A.2.1.1–A.2.1.6 del IS Code 2008 Cap.2 - gz_curve_widget.py: GZCurveWidget QPainter — curva cian, áreas sombreadas, líneas IMO, marcador AVS, tabla PASS/FAIL integrada - main_window.py: GZCurveWidget en MOD_STABILITY, _compute_and_show_gz, _on_show_stability conectado al ribbon - dark.qss: estilos GZCurveWidget - test_module3_stability.py: 33 tests S-01..S-28 (315 total, todos pasan) Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
580 lines
20 KiB
Python
580 lines
20 KiB
Python
"""
|
||
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
|