# Copyright (C) 2026 Antoine ANCEAU
#
# This file is part of navaltoolbox.
#
# navaltoolbox is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.

"""
Tests for navaltoolbox Python bindings.

These tests mirror the Rust integration tests to verify
Python bindings work correctly.
"""

import pytest
from pathlib import Path
import tempfile


class TestHullLoading:
    """Tests for Hull STL file loading."""

    def test_load_normal_stl(self):
        """Test loading a normal STL file (DTMB 5415)."""
        from navaltoolbox import Hull

        stl_path = (
            Path(__file__).parent.parent.parent
            / "rust"
            / "tests"
            / "data"
            / "dtmb5415.stl"
        )

        if not stl_path.exists():
            pytest.skip(f"DTMB 5415 STL not found at {stl_path}")

        hull = Hull(str(stl_path))
        num_triangles = hull.num_triangles()
        bounds = hull.get_bounds()

        # Verify hull loaded successfully
        assert num_triangles > 1000, (
            f"Expected > 1000 triangles, got {num_triangles}"
        )

        # Verify bounds are reasonable for DTMB 5415
        loa = bounds[1] - bounds[0]
        boa = bounds[3] - bounds[2]

        assert 140.0 < loa < 160.0, f"LOA should be ~142-153m, got {loa:.2f}m"
        assert 15.0 < boa < 25.0, f"BOA should be ~19-20m, got {boa:.2f}m"

    def test_load_box_stl(self):
        """Test loading a box STL file."""
        from navaltoolbox import Hull

        stl_path = (
            Path(__file__).parent.parent.parent
            / "rust"
            / "tests"
            / "data"
            / "box_10x10.stl"
        )

        if not stl_path.exists():
            pytest.skip(f"Box STL not found at {stl_path}")

        hull = Hull(str(stl_path))
        num_triangles = hull.num_triangles()

        # Box should have a reasonable number of triangles
        assert num_triangles > 0, "Box hull should have triangles"

    def test_empty_stl_file_error(self):
        """Test that loading an empty file raises an appropriate error."""
        from navaltoolbox import Hull

        with tempfile.NamedTemporaryFile(suffix=".stl", delete=False) as f:
            temp_path = f.name

        try:
            # Attempt to load empty file should raise an error
            with pytest.raises(Exception) as exc_info:
                _ = Hull(temp_path)

            error_msg = str(exc_info.value).lower()
            # Check that error message is meaningful
            is_meaningful = (
                "empty" in error_msg
                or "stl" in error_msg
                or "failed" in error_msg
            )
            assert is_meaningful, (
                f"Error message should mention empty/STL/failed, "
                f"got: {exc_info.value}"
            )

        finally:
            Path(temp_path).unlink(missing_ok=True)


class TestHullThickness:
    """Tests for hull thickness feature."""

    def test_thickness_properties(self):
        from navaltoolbox import Hull, Vessel
        hull = Hull.from_box(length=10.0, breadth=5.0, depth=2.0)

        # Test default is None
        assert hull.thickness is None

        # Test setter
        hull.thickness = 0.015
        assert hull.thickness == 0.015

        # Test Vessel methods
        vessel = Vessel(hull)
        assert vessel.get_hull_thickness(0) == 0.015
        vessel.set_hull_thickness(0, 0.02)
        assert vessel.get_hull_thickness(0) == 0.02

    def test_thickness_volume_addition(self):
        from navaltoolbox import Hull, Vessel, HydrostaticsCalculator
        hull = Hull.from_box(length=10.0, breadth=10.0, depth=10.0)
        vessel = Vessel(hull)
        calc1 = HydrostaticsCalculator(vessel, 1000.0)

        # Draft of 5m -> Wetted surface area = Bottom + 4 sides
        # Bottom = 10 * 10 = 100
        # 4 Sides = 4 * (10 * 5) = 200
        # Total WSA = 300 m^2
        state_no_thickness = calc1.from_draft(5.0)
        assert state_no_thickness.thickness_volume == 0.0

        vessel.set_hull_thickness(0, 0.01)
        calc2 = HydrostaticsCalculator(vessel, 1000.0)
        state_with_thickness = calc2.from_draft(5.0)

        expected_added_volume = 300.0 * 0.01
        actual_added = state_with_thickness.thickness_volume
        assert abs(actual_added - expected_added_volume) < 1e-4

        expected_total = state_no_thickness.volume + expected_added_volume
        assert abs(state_with_thickness.volume - expected_total) < 1e-4

    def test_thickness_equivalence_gz_and_hydro(self):
        """Test realistic hull thickness on displacement and GZ."""
        from navaltoolbox import (
            Hull,
            Vessel,
            HydrostaticsCalculator,
            StabilityCalculator,
        )

        # 10m x 4m x 5m box.
        hull_thick = Hull.from_box(length=10.0, breadth=4.0, depth=5.0)
        vessel_thick = Vessel(hull_thick)
        vessel_thick.set_hull_thickness(0, 0.015)  # 15mm realistic thickness

        hull_ref = Hull.from_box(length=10.0, breadth=4.0, depth=5.0)
        vessel_ref = Vessel(hull_ref)

        calc_thick = HydrostaticsCalculator(vessel_thick, 1000.0)

        # For a 10x4x5 box at 2m draft:
        # Bare Volume = 10 * 4 * 2 = 80 m^3
        # WSA = 10*4 (bottom) + 2*(10*2) + 2*(4*2) = 40 + 40 + 16 = 96 m^2
        # Added Volume = 96 * 0.015 = 1.44 m^3
        # Total Volume = 81.44 m^3
        state_thick = calc_thick.from_draft(2.0)
        assert abs(state_thick.volume - 81.44) < 1e-4

        # Test GZ equivalence at same target displacement.
        # This ensures both models submerge to roughly the same draft.
        disp = 81440.0
        cog = (5.0, 0.0, 1.0)  # center of 10m box
        heels = [0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0]

        stab_thick = StabilityCalculator(vessel_thick, 1000.0)
        stab_ref = StabilityCalculator(vessel_ref, 1000.0)

        gz_thick = stab_thick.gz_curve(disp, cog, heels)
        gz_ref = stab_ref.gz_curve(disp, cog, heels)

        # The thickness model adds volume but does not widen the waterplane.
        # For a realistic 15mm plate, the GZ curve should be nearly identical.
        for pt_t, pt_r in zip(gz_thick.points(), gz_ref.points()):
            assert abs(pt_t[3] - pt_r[3]) < 0.015, f"GZ mismatch at {pt_t[0]}°"

    def test_contact_area_wetted_surface(self):
        """Test contact surface detection between two adjacent boxes."""
        from navaltoolbox import Hull, Vessel, HydrostaticsCalculator

        # Box 1: 10x5x10, shifted to Y=-2.5
        hull1 = Hull.from_box(length=10.0, breadth=5.0, depth=10.0)
        hull1.transform((0.0, -2.5, 0.0), (0.0, 0.0, 0.0), (0.0, 0.0, 0.0))
        # Box 2: 10x5x10, shifted to Y=+2.5
        hull2 = Hull.from_box(length=10.0, breadth=5.0, depth=10.0)
        hull2.transform((0.0, 2.5, 0.0), (0.0, 0.0, 0.0), (0.0, 0.0, 0.0))

        # Together they form a 10x10x10 box
        # Face at Y=0 is shared (10x10)

        vessel = Vessel.from_hulls([hull1, hull2])
        vessel.set_hull_thickness(0, 0.1)
        vessel.set_hull_thickness(1, 0.1)

        calc = HydrostaticsCalculator(vessel, 1000.0)

        # Submerge to 5m draft
        # Each box has submerged surface = Bottom(50) + 2*Ends(50) +
        # 2*Sides(100) = 200
        # The shared side (10x5=50) should be excluded twice (once per hull)
        state = calc.from_draft(5.0)

        # Total raw wetted surface of two separate boxes = 400
        # Contact surface should be detected as 50
        assert abs(state.contact_surface_area - 50.0) < 1e-4

    def test_precomputed_contact_surfaces(self):
        """Test the pre-computation of contact surfaces between hulls."""
        from navaltoolbox import Hull, Vessel, HydrostaticsCalculator

        # Box 1: 10x5x10, shifted to Y=-2.5
        hull1 = Hull.from_box(length=10.0, breadth=5.0, depth=10.0)
        hull1.transform((0.0, -2.5, 0.0), (0.0, 0.0, 0.0), (0.0, 0.0, 0.0))
        # Box 2: 10x5x10, shifted to Y=+2.5
        hull2 = Hull.from_box(length=10.0, breadth=5.0, depth=10.0)
        hull2.transform((0.0, 2.5, 0.0), (0.0, 0.0, 0.0), (0.0, 0.0, 0.0))

        # create multi-hull vessel, which auto-computes contact surfaces
        vessel = Vessel.from_hulls([hull1, hull2])

        assert vessel.has_contact_surfaces()
        assert vessel.num_contact_surfaces() == 1

        surfaces = vessel.get_contact_surfaces()
        assert len(surfaces) == 1

        cs = surfaces[0]
        assert cs.hull_i == 0
        assert cs.hull_j == 1
        # The two hulls touch on a 10x10 surface, so 100 m^2
        assert abs(cs.total_area - 100.0) < 1e-4
        assert cs.num_faces_i > 0
        assert cs.num_faces_j > 0

        # Submerge to 5m draft and check contact surface calculation (50 m^2)
        calc = HydrostaticsCalculator(vessel, 1000.0)
        state = calc.from_draft(5.0)
        assert abs(state.contact_surface_area - 50.0) < 1e-4

        # Clear surfaces
        vessel.clear_contact_surfaces()
        assert not vessel.has_contact_surfaces()
        assert vessel.num_contact_surfaces() == 0

        # Run hydrostatics again. Fallback should be identical
        state2 = calc.from_draft(5.0)
        assert abs(state2.contact_surface_area - 50.0) < 1e-4

        # Re-compute
        vessel.compute_contact_surfaces()
        assert vessel.has_contact_surfaces()


class TestTank:
    """Tests for Tank class."""

    def test_box_tank_volume(self):
        """Test box tank volume calculation."""
        from navaltoolbox import Tank

        # 10 x 5 x 2 = 100 m³
        tank = Tank.from_box("Test", 0.0, 10.0, 0.0, 5.0, 0.0, 2.0, 1000.0)

        expected_volume = 100.0
        error = abs(tank.total_volume - expected_volume)

        assert error < 1.0, (
            f"Tank volume error: expected {expected_volume}, "
            f"got {tank.total_volume}"
        )

    def test_tank_fill_level(self):
        """Test tank fill level management."""
        from navaltoolbox import Tank

        tank = Tank.from_box("Test", 0.0, 10.0, 0.0, 5.0, 0.0, 2.0, 1000.0)

        # Set 50% fill
        tank.fill_percent = 50.0

        assert abs(tank.fill_level - 0.5) < 1e-6
        assert abs(tank.fill_percent - 50.0) < 1e-6

        # Fill volume should be ~50 m³
        assert (
            abs(tank.fill_volume - 50.0) < 1.0
        ), f"Fill volume should be ~50, got {tank.fill_volume}"

        # Fluid mass should be ~50000 kg
        assert (
            abs(tank.fluid_mass - 50000.0) < 100.0
        ), f"Fluid mass should be ~50000, got {tank.fluid_mass}"

    def test_tank_center_of_gravity(self):
        """Test tank center of gravity calculation."""
        from navaltoolbox import Tank

        tank = Tank.from_box("Test", 0.0, 10.0, -2.5, 2.5, 0.0, 2.0, 1000.0)
        tank.fill_percent = 50.0

        cog = tank.center_of_gravity

        # For a box tank at 50% fill:
        # x should be at center (5.0)
        # y should be at center (0.0)
        # z should be at half of fill height
        assert abs(cog[0] - 5.0) < 0.1, f"CoG x should be ~5.0, got {cog[0]}"
        assert abs(cog[1]) < 0.1, f"CoG y should be ~0.0, got {cog[1]}"
        assert cog[2] < 1.0, f"CoG z should be < 1.0 (half fill), got {cog[2]}"

    def test_tank_free_surface_moment(self):
        """Test free surface moment calculation."""
        from navaltoolbox import Tank

        tank = Tank.from_box("Test", 0.0, 10.0, -2.5, 2.5, 0.0, 2.0, 1000.0)
        tank.fill_percent = 50.0

        fsm_t = tank.free_surface_moment_t
        fsm_l = tank.free_surface_moment_l

        assert fsm_t > 0.0, "FSM_t should be > 0 at partial fill"
        assert fsm_l > 0.0, "FSM_l should be > 0 at partial fill"

        # For box tank:
        # I_t = L * B³ / 12 = 10 * 5³ / 12 ≈ 104.17
        # I_l = B * L³ / 12 = 5 * 10³ / 12 ≈ 416.67
        assert abs(fsm_t - 104.17) < 5.0, (
            f"FSM_t should be ~104.17, got {fsm_t}"
        )
        assert abs(fsm_l - 416.67) < 20.0, (
            f"FSM_l should be ~416.67, got {fsm_l}"
        )

    def test_tank_no_free_surface_when_full(self):
        """Test no free surface when tank is full."""
        from navaltoolbox import Tank

        tank = Tank.from_box("Test", 0.0, 10.0, -2.5, 2.5, 0.0, 2.0, 1000.0)
        tank.fill_percent = 100.0

        assert tank.free_surface_moment_t == 0.0, "FSM_t should be 0 when full"
        assert tank.free_surface_moment_l == 0.0, "FSM_l should be 0 when full"

    def test_tank_permeability(self):
        """Test tank permeability factor."""
        from navaltoolbox import Tank

        # 10 x 5 x 2 = 100 m³
        tank = Tank.from_box("Test", 0.0, 10.0, 0.0, 5.0, 0.0, 2.0, 1000.0)

        # Default permeability should be 1.0
        assert abs(tank.permeability - 1.0) < 1e-6

        # Initial state at 50% fill
        tank.fill_percent = 50.0
        orig_volume = tank.fill_volume
        orig_mass = tank.fluid_mass
        orig_fsm_t = tank.free_surface_moment_t

        # Apply 90% permeability
        tank.permeability = 0.9

        # Property should be preserved
        assert abs(tank.permeability - 0.9) < 1e-6

        # Values should be 90% of original
        assert abs(tank.fill_volume - orig_volume * 0.9) < 1e-6
        assert abs(tank.fluid_mass - orig_mass * 0.9) < 1e-6
        assert abs(tank.free_surface_moment_t - orig_fsm_t * 0.9) < 1e-6


class TestVessel:
    """Tests for Vessel class."""

    def test_vessel_from_tank_manager(self):
        """Test vessel tank management."""
        from navaltoolbox import Vessel, Tank, Hull

        # We need a hull to create a vessel - use the DTMB5415
        stl_path = (
            Path(__file__).parent.parent.parent
            / "rust"
            / "tests"
            / "data"
            / "dtmb5415.stl"
        )
        if not stl_path.exists():
            pytest.skip("DTMB5415 STL file not found")

        hull = Hull(str(stl_path))
        vessel = Vessel(hull)

        assert vessel.num_tanks() == 0

        # Add a tank
        tank = Tank.from_box("Fuel1", 50.0, 60.0, -5.0, 5.0, 0.0, 3.0, 850.0)
        vessel.add_tank(tank)

        assert vessel.num_tanks() == 1


class TestDTMB5415:
    """DTMB5415 validation tests."""

    # Reference GZ data from Ariffin (2017) PhD Thesis Figure 3.15
    REFERENCE_DATA = {
        0: 0.000,
        5: 0.171,
        10: 0.339,
        15: 0.505,
        20: 0.674,
        25: 0.848,
        30: 0.993,
        35: 1.069,
        40: 1.077,
        45: 1.025,
        50: 0.924,
        55: 0.789,
        60: 0.625,
    }

    # Loading condition
    DISPLACEMENT = 8635000.0  # kg
    LCG = 71.670
    TCG = 0.0
    VCG = 7.555

    @pytest.fixture
    def dtmb5415_vessel(self):
        """Load DTMB5415 hull and create vessel."""
        from navaltoolbox import Hull, Vessel

        stl_path = (
            Path(__file__).parent.parent.parent
            / "rust"
            / "tests"
            / "data"
            / "dtmb5415.stl"
        )
        if not stl_path.exists():
            pytest.skip("DTMB5415 STL file not found")

        hull = Hull(str(stl_path))
        return Vessel(hull)

    def test_dtmb5415_hull_loads(self):
        """Test DTMB5415 hull loads correctly."""
        from navaltoolbox import Hull

        stl_path = (
            Path(__file__).parent.parent.parent
            / "rust"
            / "tests"
            / "data"
            / "dtmb5415.stl"
        )
        if not stl_path.exists():
            pytest.skip("DTMB5415 STL file not found")

        hull = Hull(str(stl_path))

        assert hull.num_triangles() > 1000, "Hull should have many triangles"

        bounds = hull.get_bounds()
        loa = bounds[1] - bounds[0]
        boa = bounds[3] - bounds[2]

        print(f"DTMB5415: LOA={loa:.2f}m, BOA={boa:.2f}m")

        assert 150.0 < loa < 160.0, f"LOA should be ~153m, got {loa}"
        assert 18.0 < boa < 22.0, f"BOA should be ~20.5m, got {boa}"

    def test_dtmb5415_volume_at_reference_draft(self, dtmb5415_vessel):
        """Test volume at reference draft."""
        from navaltoolbox import HydrostaticsCalculator

        calc = HydrostaticsCalculator(dtmb5415_vessel, 1025.0)
        state = calc.from_draft(6.15, 0.0, 0.0, self.VCG)

        print(f"DTMB5415 at T=6.15m: Volume={state.volume:.1f}m³")

        # Reference: ~8424 m³
        assert 7500.0 < state.volume < 9000.0, (
            f"Volume should be ~8424m³, got {state.volume:.1f}"
        )

    def test_dtmb5415_gz_curve_accuracy(self, dtmb5415_vessel):
        """Test GZ curve accuracy against reference data."""
        from navaltoolbox import StabilityCalculator

        calc = StabilityCalculator(dtmb5415_vessel, 1025.0)

        heels = list(self.REFERENCE_DATA.keys())
        cog = (self.LCG, self.TCG, self.VCG)
        curve = calc.gz_curve(self.DISPLACEMENT, cog, heels)

        tolerance = 0.10  # 10cm

        passed = 0
        total = len(heels)

        print("\nDTMB5415 GZ Validation:")
        print("Heel      Calc GZ    Ref GZ     Error")
        print("-" * 45)

        for heel, gz in zip(curve.heels(), curve.values()):
            ref_gz = self.REFERENCE_DATA.get(int(heel), 0)
            error = abs(gz - ref_gz)

            status = "✓" if error < tolerance else "✗"
            if error < tolerance:
                passed += 1

            print(
                f"{heel:5.1f}°    {gz:7.3f}m   "
                f"{ref_gz:7.3f}m   {error:.3f}m {status}"
            )

        pass_rate = passed / total
        print(f"\nPassed: {passed}/{total} ({pass_rate * 100:.0f}%)")

        assert pass_rate >= 0.7, (
            f"At least 70% should pass, got {pass_rate * 100:.0f}%"
        )

    def test_dtmb5415_max_gz_location(self, dtmb5415_vessel):
        """Test maximum GZ location."""
        from navaltoolbox import StabilityCalculator

        calc = StabilityCalculator(dtmb5415_vessel, 1025.0)

        heels = list(range(0, 65, 5))
        cog = (self.LCG, self.TCG, self.VCG)
        curve = calc.gz_curve(self.DISPLACEMENT, cog, heels)

        # Find max GZ
        gz_values = curve.values()
        max_gz = max(gz_values)
        max_idx = gz_values.index(max_gz)
        max_heel = heels[max_idx]

        print(
            f"DTMB5415 Max GZ: {max_gz:.3f}m at {max_heel}° "
            f"(Ref: ~1.08m at ~38°)"
        )

        # Max should be between 30-50°
        assert 30.0 <= max_heel <= 50.0, f"Max GZ at wrong angle: {max_heel}°"

        # Max GZ should be ~1.0m
        assert 0.7 < max_gz < 1.5, f"Max GZ should be ~1.0m, got {max_gz:.3f}m"

    def test_kn_curve_binding(self, dtmb5415_vessel):
        """Test KN curve binding with multi-displacement support."""
        from navaltoolbox import StabilityCalculator

        calc = StabilityCalculator(dtmb5415_vessel, 1025.0)
        heels = [0, 10, 20]

        # Test with multiple displacements
        displacements = [self.DISPLACEMENT, self.DISPLACEMENT * 0.9]
        curves = calc.kn_curve(
            displacements, heels, lcg=self.LCG, tcg=self.TCG
        )

        assert len(curves) == 2, "Should return one curve per displacement"

        # Verify consistent results with GZ(VCG=0)
        kn_curve_0 = curves[0]

        # Compare with gz_curve(VCG=0)
        gz_curve_0 = calc.gz_curve(
            self.DISPLACEMENT, (self.LCG, self.TCG, 0.0), heels
        )

        for kn_val, gz_val in zip(kn_curve_0.values(), gz_curve_0.values()):
            assert abs(kn_val - gz_val) < 1e-6, "KN should equal GZ(VCG=0)"

    def test_fixed_trim_differs_from_free_trim(self, dtmb5415_vessel):
        """Test that fixed_trim param actively changes convergence behavior."""
        from navaltoolbox import StabilityCalculator

        calc = StabilityCalculator(dtmb5415_vessel, 1025.0)

        # We need a loading condition where LCG causes a large trim normally.
        # dtmb LCG is ~71m. Set LCG to 0.0 (the AP), ship will heavily trim.
        heels = [0.0, 10.0, 20.0]
        cog_bad_trim = (0.0, 0.0, 6.5)

        # 1. Free trim (should result in huge trim)
        curve_free = calc.gz_curve(self.DISPLACEMENT, cog_bad_trim, heels)

        # 2. Fixed trim (forced 0.0)
        curve_fixed = calc.gz_curve(
            self.DISPLACEMENT, cog_bad_trim, heels, fixed_trim=0.0
        )

        assert len(curve_free.get_stability_points()) == len(
            curve_fixed.get_stability_points()
        )

        # The points differ drastically because fixed trim ignores LCG mismatch
        point_free = curve_free.get_stability_points()[1]  # 10 degrees
        point_fixed = curve_fixed.get_stability_points()[1]

        # Free trim should be non-zero
        assert abs(point_free.trim) > 5.0, \
            "Free trim should be huge (unbalanced LCG)"

        # Fixed trim MUST remain exactly 0.0
        assert abs(point_fixed.trim) < 1e-6, \
            "Fixed trim should remain 0.0 exactly"

        # The resulting GZ should be different
        assert abs(point_free.gz - point_fixed.gz) > 0.1, \
            "GZ should differ based on trim"


class TestWallSidedFormula:
    """Wall-sided GZ formula validation (box barge)."""

    def test_wall_sided_gz(self):
        """
        Validate GZ against exact Wall-Sided Formula.

        Box: L=10, B=10, D=10.
        Draft T = 5.0m.
        KB = T/2 = 2.50m.
        BM = B² / 12T = 100 / 60 = 1.6667m.
        KM = KB + BM = 4.1667m.
        KG = 2.0m.
        GM = KM - KG = 2.1667m.

        Formula: GZ = (GM + 0.5 * BM * tan²(φ)) * sin(φ)
        """
        import math
        from navaltoolbox import Hull, Vessel, StabilityCalculator

        # Create box hull using from_box
        hull = Hull.from_box(10.0, 10.0, 10.0)
        vessel = Vessel(hull)
        stab_calc = StabilityCalculator(vessel, water_density=1025.0)

        # Parameters
        L, B = 10.0, 10.0  # D not used after box creation
        T = 5.0  # Draft
        KG = 2.0  # VCG

        # Calculated values
        KB = T / 2.0  # 2.5
        BM = (B ** 2) / (12.0 * T)  # 1.6667
        GM = KB + BM - KG  # 2.1667

        # Displacement for T=5m
        volume = L * B * T
        displacement = volume * 1025.0

        # COG: centered at (L/2, 0, KG)
        cog = (L / 2.0, 0.0, KG)

        # Calculate GZ curve
        heels = [0.0, 10.0, 20.0, 30.0]
        curve = stab_calc.gz_curve(displacement, cog, heels)

        # Validate against wall-sided formula
        # points() returns [(heel, draft, trim, gz), ...]
        for (heel, draft, trim, gz) in curve.points():
            heel_rad = math.radians(heel)
            tan_phi = math.tan(heel_rad)
            sin_phi = math.sin(heel_rad)

            # Wall-sided formula: GZ = (GM + 0.5 * BM * tan²φ) * sin(φ)
            gz_theoretical = (GM + 0.5 * BM * tan_phi ** 2) * sin_phi

            # Allow some tolerance (numerical integration vs analytical)
            tolerance = 0.05  # 5cm
            error = abs(gz - gz_theoretical)

            assert error < tolerance, (
                f"At heel={heel}°: "
                f"calculated GZ={gz:.4f}m, "
                f"theoretical={gz_theoretical:.4f}m, "
                f"error={error:.4f}m"
            )


class TestWaterplaneCalculations:
    """Waterplane property validation tests."""

    def test_dtmb5415_metacentric_height(self):
        """
        Validate GMT and KMt against DTMB 5415 reference data.

        Reference (Semarak paper, Table 3):
        - Draft: 6.15 m
        - VCG: 7.555 m
        - KMt: 9.493 m
        - GMT: 1.95 m (KMt - VCG)

        Tolerance: ±0.1 m
        """
        from navaltoolbox import Hull, Vessel, HydrostaticsCalculator

        stl_path = (
            Path(__file__).parent.parent.parent
            / "rust"
            / "tests"
            / "data"
            / "dtmb5415.stl"
        )

        if not stl_path.exists():
            pytest.skip(f"DTMB 5415 STL not found at {stl_path}")

        hull = Hull(str(stl_path))
        vessel = Vessel(hull)
        calc = HydrostaticsCalculator(vessel, 1025.0)

        # Calculate at reference draft with VCG
        state = calc.from_draft(6.15, vcg=7.555)

        # Calculate KMt
        kmt = state.vcb + state.bmt
        gmt = state.gmt

        # Reference values
        ref_kmt = 9.493
        ref_gmt = 1.95
        tolerance = 0.1

        # Validate KMt
        kmt_error = abs(kmt - ref_kmt)
        assert kmt_error < tolerance, (
            f"KMt error {kmt_error:.3f}m exceeds tolerance {tolerance}m. "
            f"Calculated: {kmt:.3f}m, Reference: {ref_kmt}m"
        )

        # Validate GMT
        gmt_error = abs(gmt - ref_gmt)
        assert gmt_error < tolerance, (
            f"GMT error {gmt_error:.3f}m exceeds tolerance {tolerance}m. "
            f"Calculated: {gmt:.3f}m, Reference: {ref_gmt}m"
        )

        print("\n✅ DTMB 5415 Waterplane Validation:")
        print(f"   KMt: {kmt:.3f}m (ref: {ref_kmt}m, error: {kmt_error:.3f}m)")
        print(f"   GMT: {gmt:.3f}m (ref: {ref_gmt}m, error: {gmt_error:.3f}m)")
        print(f"   Awp: {state.waterplane_area:.1f}m²")

    def test_waterplane_properties_exposed(self):
        """Verify waterplane properties are accessible via bindings."""
        from navaltoolbox import Hull, Vessel, HydrostaticsCalculator

        stl_path = (
            Path(__file__).parent.parent.parent
            / "rust"
            / "tests"
            / "data"
            / "dtmb5415.stl"
        )

        if not stl_path.exists():
            pytest.skip("DTMB 5415 STL not found")

        hull = Hull(str(stl_path))
        vessel = Vessel(hull)
        calc = HydrostaticsCalculator(vessel, 1025.0)

        state = calc.from_draft(6.15, vcg=7.555)

        # Check that all waterplane properties exist and are reasonable
        assert hasattr(state, "waterplane_area")
        assert hasattr(state, "lcf")
        assert hasattr(state, "bmt")
        assert hasattr(state, "bml")
        assert hasattr(state, "lwl")
        assert hasattr(state, "bwl")

        # Values should be positive and reasonable
        assert state.waterplane_area > 1000.0, "Awp should be > 1000 m²"
        assert state.bmt > 0.0, "BMt should be positive"
        assert state.bml > 0.0, "BMl should be positive"
        assert state.lwl > 100.0, "LWL should be > 100m"
        assert state.bwl > 10.0, "BWL should be > 10m"

    def test_gmt_gml_optional_without_vcg(self):
        """Verify GMT/GML are None when VCG not provided."""
        from navaltoolbox import Hull, Vessel, HydrostaticsCalculator

        stl_path = (
            Path(__file__).parent.parent.parent
            / "rust"
            / "tests"
            / "data"
            / "dtmb5415.stl"
        )

        if not stl_path.exists():
            pytest.skip("DTMB 5415 STL not found")

        hull = Hull(str(stl_path))
        vessel = Vessel(hull)
        calc = HydrostaticsCalculator(vessel, 1025.0)

        # Calculate without VCG
        state = calc.from_draft(6.15)

        # GMT/GML should be None
        assert state.gmt is None, "GMT should be None when VCG not provided"
        assert state.gml is None, "GML should be None when VCG not provided"
        assert state.cog is None, "COG should be None when VCG not provided"

        # But waterplane properties should still be calculated
        assert state.waterplane_area > 0.0
        assert state.bmt > 0.0
        assert state.bml > 0.0


class TestFreeSurfaceCorrection:
    """Free surface correction validation tests."""

    def test_calc_fsc_box_tank(self):
        """
        Validate FSC calculation for a box tank.

        Tank: 5m x 5m x 2m
        FSM_t = L * B³ / 12 = 5 * 125 / 12 = 52.083 m⁴
        Displacement (Mass) calculated from box hull.
        FSC = FSM / Displacement
        """
        from navaltoolbox import Hull, Vessel, Tank, HydrostaticsCalculator

        # 1. Setup Vessel
        stl_path = (
            Path(__file__).parent.parent.parent
            / "rust"
            / "tests"
            / "data"
            / "box_10x10.stl"
        )
        if not stl_path.exists():
            # Fallback to DTMB if box not available
            stl_path = (
                Path(__file__).parent.parent.parent
                / "rust"
                / "tests"
                / "data"
                / "dtmb5415.stl"
            )

        if not stl_path.exists():
            pytest.skip("No suitable STL found for FSC test")

        hull = Hull(str(stl_path))
        vessel = Vessel(hull)

        # 2. Add Tank
        tank = Tank.from_box(
            name="TestTank",
            x_min=0.0,
            x_max=5.0,
            y_min=-2.5,
            y_max=2.5,
            z_min=0.0,
            z_max=2.0,
            fluid_density=1000.0,
        )
        tank.fill_percent = 50.0
        vessel.add_tank(tank)

        expected_fsm = 5.0 * 5.0**3 / 12.0  # 52.083
        assert abs(tank.free_surface_moment_t - expected_fsm) < 0.1

        # 3. Calculate
        calc = HydrostaticsCalculator(vessel, 1025.0)
        state = calc.from_draft(5.0, vcg=5.0)

        # 4. Validate
        gmt_wet = state.gmt
        gmt_dry = state.gmt_dry

        # Check dry calculation
        kmt = state.vcb + state.bmt
        expected_dry = kmt - 5.0
        assert abs(gmt_dry - expected_dry) < 1e-6, "GMT dry mismatch"

        # Check wet calculation
        # FSC = (FSM * fluid_density) / displacement
        fluid_density = 1000.0  # Tank fluid density in kg/m³
        expected_fsc = (
            tank.free_surface_moment_t * fluid_density / state.displacement
        )
        actual_fsc = gmt_dry - gmt_wet

        assert abs(actual_fsc - expected_fsc) < 1e-6, (
            f"FSC mismatch: expected {expected_fsc}, got {actual_fsc}"
        )
        assert gmt_wet < gmt_dry, "Wet GMT should be less than dry GMT"


class TestFreeboardImplementation:
    """Freeboard property exposed correctly in StabilityPoint."""

    def test_freeboard_with_deck_edge(self):
        from navaltoolbox import (
            Hull,
            Vessel,
            DeckEdge,
            DeckEdgeSide,
            StabilityCalculator,
        )

        vessel = Vessel(Hull.from_box(10.0, 10.0, 10.0))
        # Add a deck edge at z=10.0, matching the hull bounds: x from 0 to 10
        edge = DeckEdge.from_points(
            "MainDeck",
            [
                (10.0, 5.0, 10.0),
                (0.0, 5.0, 10.0),
                (0.0, -5.0, 10.0),
                (10.0, -5.0, 10.0)
            ],
            DeckEdgeSide.both()
        )
        vessel.add_deck_edge(edge)

        calc = StabilityCalculator(vessel, 1025.0)
        # Displacement for draft = 5.0m is 10 * 10 * 5.0 * 1025.0 = 512500.0 kg
        heels = [0.0, 10.0, 20.0]
        # Set LCG = 5.0 to align exactly with centroid and maintain 0 trim
        curve = calc.gz_curve(512500.0, (5.0, 0.0, 5.0), heels)

        points = curve.get_stability_points()
        for point in points:
            assert point.freeboard is not None

        # At 0 heel and 0 trim, freeboard from Deck 10 should
        # be exactly 10.0 - draft
        expected_freeboard = 10.0 - points[0].draft
        assert abs(points[0].freeboard - expected_freeboard) < 1e-4

    def test_freeboard_without_deck_edge(self):
        from navaltoolbox import Hull, Vessel, StabilityCalculator

        vessel = Vessel(Hull.from_box(10.0, 10.0, 10.0))
        calc = StabilityCalculator(vessel, 1025.0)
        heels = [0.0]
        curve = calc.gz_curve(512500.0, (0.0, 0.0, 5.0), heels)

        points = curve.get_stability_points()
        assert points[0].freeboard is None


if __name__ == "__main__":
    pytest.main([__file__, "-v"])