Stability Tutorial
==================

This tutorial covers ship stability calculations, including GZ curves
(righting lever curves) and KN curves (cross curves of stability).

Introduction
------------

Ship stability analysis determines how a vessel responds to heeling forces
and whether it will return to an upright position.

Key concepts:

* **GZ (Righting Lever)**: The horizontal distance between G and B when heeled
* **GZ Curve**: Plot of GZ vs heel angle - shows stability characteristics
* **Range of Stability**: Maximum heel angle with positive righting moment
* **KN Curves**: Cross curves assuming KG = 0, used for multiple loading conditions

The Righting Lever (GZ)
-----------------------

When a ship heels, the center of buoyancy (B) moves. The righting lever (GZ) is:

.. math::

   GZ = KN - KG \\cdot \\sin(\\phi)

Where:

* KN = righting lever assuming KG = 0
* KG = height of center of gravity above keel
* φ = heel angle

Setting Up
----------

First, load your hull and create the stability calculator:

.. code-block:: python

   from pynavaltoolbox import Hull, Vessel
   from pynavaltoolbox.stability import StabilityCalculator

   hull_geom = Hull("path/to/hull.stl")
   vessel = Vessel(hull_geom)
   
   stability = StabilityCalculator(vessel, water_density=1025.0)

Calculating a GZ Curve
----------------------

The GZ curve is calculated for a specific displacement and center of gravity:

.. code-block:: python

   # Define loading condition
   displacement = 5000000.0  # kg (5000 tonnes)
   cog = (70.0, 0.0, 10.0)  # LCG, TCG, VCG in meters

   # Calculate GZ curve from 0° to 90° in 5° increments
   gz_curve = stability.calculate_gz_curve(
       displacement_mass=displacement,
       cog=cog,
       heels=list(range(0, 91, 5))
   )

   # Print results
   print("Heel(°)\tGZ(m)\tDraft(m)\tTrim(°)")
   print("-" * 50)
   for point in gz_curve.points:
       print(f"{point.heel}\t{point.gz:.3f}\t{point.draft:.2f}\t\t{point.trim:.2f}")

Analyzing the GZ Curve
----------------------

The :class:`~pynavaltoolbox.stability.StabilityCurve` provides methods for analysis:

.. code-block:: python

   # Key stability characteristics
   max_gz = gz_curve.get_max_gz()
   angle_max_gz = gz_curve.get_angle_of_max_gz()
   range_stability = gz_curve.get_range_of_stability()

   print(f"Maximum GZ: {max_gz:.3f} m")
   print(f"Angle of Maximum GZ: {angle_max_gz:.1f}°")
   print(f"Range of Stability: {range_stability:.1f}°")

   # Access raw data for plotting
   heels = gz_curve.heels  # numpy array
   gz_values = gz_curve.gz_values  # numpy array

Plotting the GZ Curve
---------------------

Create a publication-quality GZ curve plot:

.. code-block:: python

   import matplotlib.pyplot as plt

   fig, ax = plt.subplots(figsize=(10, 6))

   ax.plot(gz_curve.heels, gz_curve.gz_values, 'b-', linewidth=2)
   ax.axhline(y=0, color='k', linestyle='-', linewidth=0.5)

   # Mark key points
   ax.axvline(x=angle_max_gz, color='r', linestyle='--', alpha=0.5)
   ax.scatter([angle_max_gz], [max_gz], color='r', s=100, zorder=5)
   ax.annotate(f'Max GZ = {max_gz:.3f} m at {angle_max_gz:.0f}°',
               xy=(angle_max_gz, max_gz),
               xytext=(angle_max_gz + 10, max_gz),
               fontsize=10)

   ax.set_xlabel('Heel Angle (degrees)', fontsize=12)
   ax.set_ylabel('GZ (m)', fontsize=12)
   ax.set_title('Righting Lever Curve', fontsize=14)
   ax.grid(True, alpha=0.3)
   ax.set_xlim(0, 90)

   plt.tight_layout()
   plt.savefig('gz_curve.png', dpi=300)
   plt.show()

Calculating KN Curves
---------------------

KN curves (cross curves of stability) are calculated assuming KG = 0:

.. code-block:: python

   # Calculate KN curve at specific displacement
   kn_curve = stability.calculate_kn_curve(
       displacement_mass=5000000.0,
       heels=[0, 5, 10, 15, 20, 30, 45, 60, 75, 90]
   )

   print("Heel(°)\tKN(m)")
   print("-" * 20)
   for point in kn_curve.points:
       print(f"{point.heel}\t{point.kn:.3f}")

Using KN Curves for Different KG Values
---------------------------------------

KN curves allow quick GZ calculation for different loading conditions:

.. code-block:: python

   import numpy as np

   # Different KG values to evaluate
   kg_values = [8.0, 10.0, 12.0]

   print("Heel\t", end="")
   for kg in kg_values:
       print(f"GZ(KG={kg})\t", end="")
   print()

   for point in kn_curve.points:
       print(f"{point.heel}°\t", end="")
       for kg in kg_values:
           gz = point.kn - kg * np.sin(np.radians(point.heel))
           print(f"{gz:.3f}\t\t", end="")
       print()

Effect of VCG on Stability
--------------------------

This example shows how moving the center of gravity affects stability:

.. code-block:: python

   import matplotlib.pyplot as plt

   fig, ax = plt.subplots(figsize=(10, 6))

   vcg_values = [8.0, 10.0, 12.0, 14.0]
   colors = ['green', 'blue', 'orange', 'red']

   for vcg, color in zip(vcg_values, colors):
       curve = stability.calculate_gz_curve(
           displacement_mass=5000000.0,
           cog=(70.0, 0.0, vcg),
           heels=list(range(0, 91, 5))
       )
       ax.plot(curve.heels, curve.gz_values, color=color,
               linewidth=2, label=f'VCG = {vcg} m')

   ax.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
   ax.set_xlabel('Heel Angle (degrees)', fontsize=12)
   ax.set_ylabel('GZ (m)', fontsize=12)
   ax.set_title('Effect of VCG on Stability', fontsize=14)
   ax.legend()
   ax.grid(True, alpha=0.3)
   ax.set_xlim(0, 90)

   plt.tight_layout()
   plt.savefig('vcg_comparison.png', dpi=300)
   plt.show()

Stability Criteria
------------------

Common stability criteria (IMO A.749) to check:

1. **Area under GZ curve** from 0° to 30°: ≥ 0.055 m·rad
2. **Area under GZ curve** from 0° to 40°: ≥ 0.090 m·rad
3. **Area under GZ curve** between 30° and 40°: ≥ 0.030 m·rad
4. **GZ at 30°**: ≥ 0.20 m
5. **Angle of maximum GZ**: ≥ 25°
6. **Initial GM**: ≥ 0.15 m

.. code-block:: python

   import numpy as np
   from scipy.integrate import trapezoid

   # Calculate area under curve (convert to radians for integration)
   heels_rad = np.radians(gz_curve.heels)
   gz = gz_curve.gz_values

   # Find indices for key angles
   idx_30 = np.argmin(np.abs(gz_curve.heels - 30))
   idx_40 = np.argmin(np.abs(gz_curve.heels - 40))

   area_0_30 = trapezoid(gz[:idx_30+1], heels_rad[:idx_30+1])
   area_0_40 = trapezoid(gz[:idx_40+1], heels_rad[:idx_40+1])
   area_30_40 = trapezoid(gz[idx_30:idx_40+1], heels_rad[idx_30:idx_40+1])

   print("=== Stability Criteria Check ===")
   print(f"Area 0-30°: {area_0_30:.4f} m·rad (min: 0.055)")
   print(f"Area 0-40°: {area_0_40:.4f} m·rad (min: 0.090)")
   print(f"Area 30-40°: {area_30_40:.4f} m·rad (min: 0.030)")
   print(f"GZ at 30°: {gz[idx_30]:.3f} m (min: 0.20)")
   print(f"Angle of max GZ: {gz_curve.get_angle_of_max_gz():.1f}° (min: 25°)")

Next Steps
----------

* See :doc:`hydrostatics` for hydrostatic calculations
* See :doc:`../api/stability` for complete API reference