feat: scipy 5-DOF optimizer, depowering, and input validation

runVPP.py

@@ -1,4 +1,4 @@
-#!/opt/miniconda3/bin/python
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 import numpy as np
 
@@ -36,7 +36,7 @@
 vpp = VPP(Yacht=YD41)
 
 vpp.set_analysis(
-    tws_range=np.arange(4.0, 22.0, 2.0), twa_range=np.linspace(30.0, 180.0, 31)
+    tws_range=np.arange(4.0, 22.0, 2.0), twa_range=np.linspace(28.0, 180.0, 39)
 )
 
 vpp.run(verbose=False)

src/AeroMod.py

@@ -7,18 +7,46 @@
 __version__ = "1.0.1"
 __email__ = "M.Lauber@soton.ac.uk"
 
-import numpy as np
-from scipy.interpolate import interp1d
-from scipy.optimize import fsolve
-from scipy.optimize import root
+import functools
+
 import matplotlib.pyplot as plt
+import numpy as np
+
 from src.UtilsMod import build_interp_func
 
 
+@functools.lru_cache(maxsize=512)
+def wind_triangle(tws, twa, vb):
+    """Analytical wind triangle: TWS/TWA/VB → (AWA, AWS) in degrees.
+
+    Uses the law of cosines on the velocity triangle formed by the true
+    wind, boat speed, and apparent wind vectors.
+    """
+    twa_rad = np.radians(twa)
+    aws = np.sqrt(tws**2 + vb**2 + 2 * tws * vb * np.cos(twa_rad))
+    if aws < 1e-12:
+        return twa, 0.0
+    cos_awa = np.clip((tws * np.cos(twa_rad) + vb) / aws, -1.0, 1.0)
+    awa = np.degrees(np.arccos(cos_awa))
+    return awa, aws
+
+
 class AeroMod(object):
     def __init__(self, Yacht, rho=1.225, mu=0.0000181):
         """
-        Initializes an Aero Model, given a set of sails
+        Aerodynamic force model.
+
+        Computes sail drive force, side force, and heeling moment from the
+        yacht's sail plan using ORC aerodynamic coefficients.
+
+        Parameters
+        ----------
+        Yacht : Yacht
+            Yacht object containing sail definitions and hull geometry.
+        rho : float, optional
+            Air density (kg/m^3). Default is 1.225 (ISA sea level).
+        mu : float, optional
+            Dynamic viscosity of air (Pa.s). Default is 1.81e-5.
         """
         # physical params
         self.rho = rho
@@ -74,7 +102,31 @@ def _measure_sails(self):
     # prototype top function in hydro mod
     def update(self, vb, phi, tws, twa, flat, RED):
         """
-        Update the aero model for current iter
+        Update aerodynamic forces for current sailing state.
+
+        Solves the wind triangle, computes sail coefficients, and projects
+        forces into the boat reference frame.
+
+        Parameters
+        ----------
+        vb : float
+            Boat speed (m/s).
+        phi : float
+            Heel angle (degrees).
+        tws : float
+            True wind speed (m/s).
+        twa : float
+            True wind angle (degrees).
+        flat : float
+            Sail flattening factor (0.62 to 1.0). Reduces lift and drag.
+        RED : float
+            Reef/reduction factor. Values > 1 apply jib furling (ftj = RED - 1),
+            values <= 1 apply mainsail reefing (rfm = RED).
+
+        Returns
+        -------
+        tuple of float
+            (Fx, Fy, Mx) — drive force (N), side force (N), heeling moment (N.m).
         """
         self.vb = max(0, vb)
         self.phi = max(0, phi)
@@ -112,7 +164,7 @@ def _compute_forces(self):
 
         # side-force is horizontal component of Fh
         self.Fy *= np.cos(np.radians(self.phi))
-        
+
         # heeling moment
         self.Mx = self.Fy * self._vce()
 
@@ -136,12 +188,15 @@ def _get_coeffs(self):
         self.cl = 0.0
         self.cd = 0.0
         kpp = 0.0
+        sail_cd = {}
 
         for sail in self.sails:
-
-            self.cl += sail.cl(self.awa) * sail.area * sail.bk
-            self.cd += sail.cd(self.awa) * sail.area * sail.bk
-            kpp += sail.cl(self.awa) ** 2 * sail.area * sail.bk * sail.kp
+            cl_i = sail.cl(self.awa)
+            cd_i = sail.cd(self.awa)
+            sail_cd[id(sail)] = cd_i
+            self.cl += cl_i * sail.area * sail.bk
+            self.cd += cd_i * sail.area * sail.bk
+            kpp += cl_i ** 2 * sail.area * sail.bk * sail.kp
 
         self.cl /= self.area
         self.cd /= self.area
@@ -156,7 +211,7 @@ def _get_coeffs(self):
         for sail in self.sails:
             if sail.type == "jib":
                 self.fcdj = (
-                    sail.bk * sail.cd(self.awa) * sail.area / (self.cd * self.area)
+                    sail.bk * sail_cd[id(sail)] * sail.area / (self.cd * self.area)
                 )
 
         # final lift and drag
@@ -170,17 +225,7 @@ def _update_windTriangle(self):
         """
         find AWS and AWA for a given TWS, TWA and VB
         """
-        _awa_ = lambda awa: self.vb * np.sin(awa / 180.0 * np.pi) - self.tws * np.sin(
-            (self.twa - awa) / 180.0 * np.pi
-        )
-        self.awa = fsolve(_awa_, self.twa)[0]
-        self.aws = np.sqrt(
-            (self.tws * np.sin(self.twa / 180.0 * np.pi)) ** 2
-            + (self.tws * np.cos(self.twa / 180.0 * np.pi) + self.vb) ** 2
-        )
-        # self.awa = np.arccos((self.tws*np.cos(np.radians(self.twa)) + self.vb) / np.sqrt((self.tws**2) + (self.vb**2) + 
-        #              2*self.tws*self.vb * np.cos(np.radians(self.twa))))
-        # self.aws = (self.tws * np.sin(np.radians(self.twa))) / np.sin(self.awa)
+        self.awa, self.aws = wind_triangle(self.tws, self.twa, self.vb)
 
 
     def _area(self):
@@ -194,7 +239,7 @@ def _area(self):
 
     def _vce(self):
         """
-        Vectical centre of effort lift/drag weigted
+        Vertical centre of effort, lift/drag weighted.
         """
         sum = 0.0
         for sail in self.sails:

src/HydroMod.py

@@ -7,14 +7,33 @@
 __version__ = "1.0.1"
 __email__ = "M.Lauber@soton.ac.uk"
 
-import numpy as np
-from scipy.interpolate import RegularGridInterpolator
 import warnings
+
 import matplotlib.pyplot as plt
+import numpy as np
+from scipy.interpolate import RegularGridInterpolator
 
 
 class HydroMod(object):
     def __init__(self, Yacht, rho=1025.0, mu=0.00119, g=9.81):
+        """
+        Hydrodynamic resistance and righting moment model.
+
+        Computes total resistance (viscous + residuary + induced), side force
+        from appendage lift, and righting moment using ORC methods and
+        interpolated resistance surfaces.
+
+        Parameters
+        ----------
+        Yacht : Yacht
+            Yacht object containing hull geometry and appendage definitions.
+        rho : float, optional
+            Seawater density (kg/m^3). Default is 1025.0.
+        mu : float, optional
+            Dynamic viscosity of seawater (Pa.s). Default is 1.19e-3.
+        g : float, optional
+            Gravitational acceleration (m/s^2). Default is 9.81.
+        """
 
         # physical parameters
         self.rho = rho
@@ -111,7 +130,7 @@ def _get_Ri(self):
             self.Teff = np.hstack((self.Teff, appendage.teff))
 
         self.Ksfj = (
-            0.5 * self.rho * self.vb ** 2 * self.cla * self.leeway / 180.0 * np.pi
+            0.5 * self.rho * self.vb ** 2 * self.cla * np.radians(self.leeway)
         )
         self.Ksf = np.sum(self.Ksfj)
 
@@ -122,24 +141,99 @@ def _get_Ri(self):
 
     def _cf(self, L):
         """
-        Flate plate turbulent boudnary layer friction coefficient.
-        Take a length scale, such that it can be used for appendags as well
+        Flat plate turbulent boundary layer friction coefficient (ITTC 1957)
+        with ITTC 1978 roughness allowance. Takes a length scale so it can
+        be used for hull and appendages.
         """
         Re = max(
             1e4, self.vb * L / self.nu
         )  # prevents dividing by zero, lowest for turbulence on plate
-        return 0.066 * (np.log10(Re) - 2.03) ** (-2)
+        cf = 0.066 * (np.log10(Re) - 2.03) ** (-2)
+        ks = self.yacht.roughness
+        if ks > 0:
+            cf += (105.0 * (ks / self.l) ** (1.0 / 3.0) - 0.64) * 1e-3
+        return cf
+
+    def _added_resistance_waves(self, twa):
+        """Added resistance in waves (simplified Gerritsma scaling).
+
+        Uses a Bretschneider-type spectrum with hull-form scaling to
+        estimate the mean added resistance from ocean waves.
+
+        Parameters
+        ----------
+        twa : float
+            True wind angle (degrees), used as wave encounter angle
+            unless ``yacht.wave_direction`` overrides it.
+
+        Returns
+        -------
+        float
+            Added resistance in Newtons. Returns 0 when Hs = 0.
+        """
+        Hs = self.yacht.Hs
+        Ts = self.yacht.Ts
+        if Hs <= 0 or Ts <= 0:
+            return 0.0
 
-    def update(self, vb, phi, leeway):
+        # wave encounter angle
+        if self.yacht.wave_direction is not None:
+            mu = np.radians(twa - self.yacht.wave_direction)
+        else:
+            mu = np.radians(twa)
+
+        # heading correction — full drag head-on, near-zero following
+        cos2_mu = np.cos(mu) ** 2
+
+        # hull-form coefficient (empirical, typical displacement hull)
+        C_aw = 6.0
+
+        # wave encounter frequency
+        omega_0 = 2.0 * np.pi / Ts
+        omega_e = abs(omega_0 - omega_0 ** 2 * self.vb * np.cos(mu) / self.g)
+
+        # resonance tuning factor (peak when encounter ≈ hull natural freq)
+        omega_n = np.sqrt(self.g / self.l)
+        r = omega_e / omega_n if omega_n > 0 else 0.0
+        f_omega = r ** 2 * np.exp(1.0 - r ** 2) if r > 0 else 0.0
+
+        Raw = C_aw * (Hs ** 2 / self.l) * (self.bwl ** 2 / self.tc) * f_omega * cos2_mu
+        # convert to Newtons (rho * g scaling)
+        Raw *= self.rho * self.g / 1000.0
+
+        return max(0.0, Raw)
+
+    def update(self, vb, phi, leeway, twa=0.0):
+        """
+        Update hydrodynamic forces for current sailing state.
+
+        Parameters
+        ----------
+        vb : float
+            Boat speed (m/s).
+        phi : float
+            Heel angle (degrees).
+        leeway : float
+            Leeway angle (degrees).
+        twa : float, optional
+            True wind angle (degrees). Used for wave encounter angle.
+
+        Returns
+        -------
+        tuple of float
+            (Fx, Fy, Mx) — total resistance (N), side force (N),
+            total righting moment (N.m).
+        """
 
         self.vb = max(0, vb)
         self.phi = max(0, phi)
-        self.leeway = max(0, leeway)
+        self.leeway = leeway
         self.lsm, self.lvr, self.btr = self.yacht.measureLSM()
         self.fn = self.vb / (np.sqrt(self.g * self.lsm))
 
         # resistance
         self.Fx = self._get_Rr() + self._get_Rv() + self._get_Ri()
+        self.Fx += self._added_resistance_waves(twa)
 
         # keel side force, calculated when _get_Ri() is called
         self.Fy = self.Ksf * np.cos(self.phi / 180.0 * np.pi)