navaltoolbox 0.9.2

// 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/>.

//! Silhouette struct for wind heeling moment calculations.

use super::loader::{
    load_csv_silhouette, load_dxf_silhouette, load_vtk_silhouette, SilhouetteLoadError,
};
use std::path::Path;

/// A 2D silhouette profile in the X-Z plane (ship lateral view).
///
/// Used for calculating wind heeling moments per IMO 2008 IS Code.
#[derive(Clone, Debug)]
pub struct Silhouette {
    /// Points defining the silhouette contour [x, y, z].
    /// Y is always 0 (X-Z plane). Stored as 3D for future extensibility.
    points: Vec<[f64; 3]>,
    /// Name of the silhouette (from filename or user-defined).
    name: String,
}

impl Silhouette {
    /// Create a new silhouette from a list of 3D points.
    pub fn new(points: Vec<[f64; 3]>, name: String) -> Self {
        let s = Self { points, name };
        s.validate();
        s
    }

    /// Validate the silhouette geometry and log warnings if issues found.
    pub fn validate(&self) {
        if self.points.len() < 3 {
            log::warn!("Silhouette '{}' has fewer than 3 points.", self.name);
        }
        if !self.is_closed() {
            log::warn!("Silhouette '{}' is not closed.", self.name);
        }
        if self.get_area() < 1e-4 {
            log::warn!(
                "Silhouette '{}' has extremely small or zero area ({:.6} m²). Check loading units or orientation.",
                self.name,
                self.get_area()
            );
        }
    }

    /// Load a silhouette from a DXF file.
    ///
    /// Extracts the first LWPOLYLINE or POLYLINE entity.
    /// If Y coordinates are non-zero, they are set to 0 with a warning.
    pub fn from_dxf(path: &Path) -> Result<Self, SilhouetteLoadError> {
        let (points, name) = load_dxf_silhouette(path)?;
        let s = Self { points, name };
        s.validate();
        Ok(s)
    }

    /// Load a silhouette from a VTK file (.vtk or .vtp).
    ///
    /// Extracts the first polyline from the PolyData.
    /// If Y coordinates are non-zero, they are set to 0 with a warning.
    pub fn from_vtk(path: &Path) -> Result<Self, SilhouetteLoadError> {
        let (points, name) = load_vtk_silhouette(path)?;
        let s = Self { points, name };
        s.validate();
        Ok(s)
    }

    /// Load a silhouette from a CSV or TXT file.
    pub fn from_csv(path: &Path) -> Result<Self, SilhouetteLoadError> {
        let (points, name) = load_csv_silhouette(path)?;
        let s = Self { points, name };
        s.validate();
        Ok(s)
    }

    /// Load a silhouette from a file (DXF or VTK) based on extension.
    pub fn from_file(path: &Path) -> Result<Self, SilhouetteLoadError> {
        let ext = path
            .extension()
            .and_then(|e| e.to_str())
            .map(|e| e.to_lowercase())
            .unwrap_or_default();

        match ext.as_str() {
            "dxf" => Self::from_dxf(path),
            "vtk" | "vtp" | "vtu" => Self::from_vtk(path),
            "csv" | "txt" => Self::from_csv(path),
            _ => Err(SilhouetteLoadError::UnsupportedFormat),
        }
    }

    /// Get the silhouette name.
    pub fn name(&self) -> &str {
        &self.name
    }

    /// Set the silhouette name.
    pub fn set_name(&mut self, name: &str) {
        self.name = name.to_string();
    }

    /// Get the points defining the silhouette.
    pub fn points(&self) -> &[[f64; 3]] {
        &self.points
    }

    /// Get number of points.
    pub fn num_points(&self) -> usize {
        self.points.len()
    }

    /// Check if the contour is closed (first point == last point).
    pub fn is_closed(&self) -> bool {
        if self.points.len() < 2 {
            return false;
        }
        let first = self.points.first().unwrap();
        let last = self.points.last().unwrap();
        let eps = 1e-6;
        (first[0] - last[0]).abs() < eps
            && (first[1] - last[1]).abs() < eps
            && (first[2] - last[2]).abs() < eps
    }

    /// Calculate the total lateral area of the silhouette (m²).
    ///
    /// Uses the shoelace formula in the X-Z plane.
    pub fn get_area(&self) -> f64 {
        if self.points.len() < 3 {
            return 0.0;
        }

        let mut area = 0.0;
        let n = self.points.len();

        for i in 0..n {
            let j = (i + 1) % n;
            let xi = self.points[i][0];
            let zi = self.points[i][2];
            let xj = self.points[j][0];
            let zj = self.points[j][2];
            area += xi * zj - xj * zi;
        }

        (area / 2.0).abs()
    }

    /// Calculate the centroid of the silhouette in the X-Z plane.
    ///
    /// Returns [x_center, z_center].
    pub fn get_centroid(&self) -> [f64; 2] {
        if self.points.len() < 3 {
            return [0.0, 0.0];
        }

        // Compute signed area (positive = CCW, negative = CW)
        let mut signed_area = 0.0;
        let n = self.points.len();
        for i in 0..n {
            let j = (i + 1) % n;
            let xi = self.points[i][0];
            let zi = self.points[i][2];
            let xj = self.points[j][0];
            let zj = self.points[j][2];
            signed_area += xi * zj - xj * zi;
        }
        signed_area /= 2.0;

        if signed_area.abs() < 1e-9 {
            return [0.0, 0.0];
        }

        let mut cx = 0.0;
        let mut cz = 0.0;
        for i in 0..n {
            let j = (i + 1) % n;
            let xi = self.points[i][0];
            let zi = self.points[i][2];
            let xj = self.points[j][0];
            let zj = self.points[j][2];
            let cross = xi * zj - xj * zi;
            cx += (xi + xj) * cross;
            cz += (zi + zj) * cross;
        }

        // Divide by signed area (not abs) so sign is handled consistently
        let factor = 1.0 / (6.0 * signed_area);
        [cx * factor, cz * factor]
    }

    /// Get the bounding box of the silhouette.
    ///
    /// Returns (x_min, x_max, z_min, z_max).
    pub fn get_bounds(&self) -> (f64, f64, f64, f64) {
        if self.points.is_empty() {
            return (0.0, 0.0, 0.0, 0.0);
        }

        let mut x_min = f64::MAX;
        let mut x_max = f64::MIN;
        let mut z_min = f64::MAX;
        let mut z_max = f64::MIN;

        for p in &self.points {
            x_min = x_min.min(p[0]);
            x_max = x_max.max(p[0]);
            z_min = z_min.min(p[2]);
            z_max = z_max.max(p[2]);
        }

        (x_min, x_max, z_min, z_max)
    }

    /// Calculate the emerged area (above waterline) in m².
    ///
    /// Clips the silhouette at the given waterline Z and returns
    /// the area of the portion above the waterline.
    pub fn get_emerged_area(&self, waterline_z: f64) -> f64 {
        let clipped = self.clip_above(waterline_z);
        Self::polygon_area(&clipped)
    }

    /// Calculate the centroid of the emerged area.
    ///
    /// Returns [x_center, z_center] of the area above waterline.
    pub fn get_emerged_centroid(&self, waterline_z: f64) -> [f64; 2] {
        let clipped = self.clip_above(waterline_z);
        Self::polygon_centroid(&clipped)
    }

    /// Calculate the submerged area (below waterline) in m².
    ///
    /// Clips the silhouette at the given waterline Z and returns
    /// the area of the portion below the waterline.
    pub fn get_submerged_area(&self, waterline_z: f64) -> f64 {
        let clipped = self.clip_below(waterline_z);
        Self::polygon_area(&clipped)
    }

    /// Calculate the centroid of the submerged area.
    ///
    /// Returns [x_center, z_center] of the area below waterline.
    /// Used for the exact Z computation per IMO 2008 IS Code §2.3.2.
    pub fn get_submerged_centroid(&self, waterline_z: f64) -> [f64; 2] {
        let clipped = self.clip_below(waterline_z);
        Self::polygon_centroid(&clipped)
    }

    /// Clip the silhouette to keep only points above the waterline.
    fn clip_above(&self, waterline_z: f64) -> Vec<[f64; 2]> {
        if self.points.len() < 2 {
            return Vec::new();
        }

        let mut result: Vec<[f64; 2]> = Vec::new();
        let n = self.points.len();

        for i in 0..n {
            let j = (i + 1) % n;
            let p1 = [self.points[i][0], self.points[i][2]];
            let p2 = [self.points[j][0], self.points[j][2]];

            let z1 = p1[1];
            let z2 = p2[1];

            // Both above
            if z1 >= waterline_z && z2 >= waterline_z {
                result.push(p1);
            }
            // p1 above, p2 below - add p1 and intersection
            else if z1 >= waterline_z && z2 < waterline_z {
                result.push(p1);
                let t = (waterline_z - z1) / (z2 - z1);
                let x_int = p1[0] + t * (p2[0] - p1[0]);
                result.push([x_int, waterline_z]);
            }
            // p1 below, p2 above - add intersection
            else if z1 < waterline_z && z2 >= waterline_z {
                let t = (waterline_z - z1) / (z2 - z1);
                let x_int = p1[0] + t * (p2[0] - p1[0]);
                result.push([x_int, waterline_z]);
            }
            // Both below - skip
        }

        result
    }

    /// Clip the silhouette to keep only points below the waterline.
    fn clip_below(&self, waterline_z: f64) -> Vec<[f64; 2]> {
        if self.points.len() < 2 {
            return Vec::new();
        }

        let mut result: Vec<[f64; 2]> = Vec::new();
        let n = self.points.len();

        for i in 0..n {
            let j = (i + 1) % n;
            let p1 = [self.points[i][0], self.points[i][2]];
            let p2 = [self.points[j][0], self.points[j][2]];

            let z1 = p1[1];
            let z2 = p2[1];

            // Both below
            if z1 <= waterline_z && z2 <= waterline_z {
                result.push(p1);
            }
            // p1 below, p2 above - add p1 and intersection
            else if z1 <= waterline_z && z2 > waterline_z {
                result.push(p1);
                let t = (waterline_z - z1) / (z2 - z1);
                let x_int = p1[0] + t * (p2[0] - p1[0]);
                result.push([x_int, waterline_z]);
            }
            // p1 above, p2 below - add intersection
            else if z1 > waterline_z && z2 <= waterline_z {
                let t = (waterline_z - z1) / (z2 - z1);
                let x_int = p1[0] + t * (p2[0] - p1[0]);
                result.push([x_int, waterline_z]);
            }
            // Both above - skip
        }

        result
    }

    /// Calculate area of a 2D polygon using shoelace formula.
    fn polygon_area(points: &[[f64; 2]]) -> f64 {
        if points.len() < 3 {
            return 0.0;
        }

        let mut area = 0.0;
        let n = points.len();

        for i in 0..n {
            let j = (i + 1) % n;
            area += points[i][0] * points[j][1] - points[j][0] * points[i][1];
        }

        (area / 2.0).abs()
    }

    /// Calculate centroid of a 2D polygon.
    ///
    /// Uses the signed-area form of the shoelace centroid formula.
    /// Handles both CW and CCW winding correctly.
    fn polygon_centroid(points: &[[f64; 2]]) -> [f64; 2] {
        if points.len() < 3 {
            return [0.0, 0.0];
        }

        let n = points.len();

        // Compute signed area (positive = CCW, negative = CW)
        let mut signed_area = 0.0;
        for i in 0..n {
            let j = (i + 1) % n;
            signed_area += points[i][0] * points[j][1] - points[j][0] * points[i][1];
        }
        signed_area /= 2.0;

        if signed_area.abs() < 1e-9 {
            return [0.0, 0.0];
        }

        let mut cx = 0.0;
        let mut cy = 0.0;
        for i in 0..n {
            let j = (i + 1) % n;
            let cross = points[i][0] * points[j][1] - points[j][0] * points[i][1];
            cx += (points[i][0] + points[j][0]) * cross;
            cy += (points[i][1] + points[j][1]) * cross;
        }

        // Divide by signed area (not abs) so sign is handled consistently
        let factor = 1.0 / (6.0 * signed_area);
        [cx * factor, cy * factor]
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn create_rectangle(x_min: f64, x_max: f64, z_min: f64, z_max: f64) -> Silhouette {
        let points = vec![
            [x_min, 0.0, z_min],
            [x_max, 0.0, z_min],
            [x_max, 0.0, z_max],
            [x_min, 0.0, z_max],
            [x_min, 0.0, z_min], // Close the polygon
        ];
        Silhouette::new(points, "rectangle".to_string())
    }

    #[test]
    fn test_rectangle_is_closed() {
        let rect = create_rectangle(0.0, 100.0, 0.0, 20.0);
        assert!(rect.is_closed());
    }

    #[test]
    fn test_rectangle_area() {
        let rect = create_rectangle(0.0, 100.0, 0.0, 20.0);
        let area = rect.get_area();
        assert!(
            (area - 2000.0).abs() < 1.0,
            "Area should be 2000 m², got {}",
            area
        );
    }

    #[test]
    fn test_emerged_area_full() {
        let rect = create_rectangle(0.0, 100.0, 0.0, 20.0);
        let emerged = rect.get_emerged_area(0.0);
        assert!(
            (emerged - 2000.0).abs() < 1.0,
            "Full emerged area should be 2000 m²"
        );
    }

    #[test]
    fn test_emerged_area_half() {
        let rect = create_rectangle(0.0, 100.0, 0.0, 20.0);
        let emerged = rect.get_emerged_area(10.0);
        assert!(
            (emerged - 1000.0).abs() < 10.0,
            "Half emerged area should be ~1000 m², got {}",
            emerged
        );
    }

    #[test]
    fn test_emerged_area_none() {
        let rect = create_rectangle(0.0, 100.0, 0.0, 20.0);
        let emerged = rect.get_emerged_area(25.0);
        assert!(emerged < 1.0, "No emerged area above z=25, got {}", emerged);
    }

    #[test]
    fn test_open_silhouette_detection() {
        let points = vec![
            [0.0, 0.0, 0.0],
            [10.0, 0.0, 0.0],
            [10.0, 0.0, 5.0],
            // Missing closure back to 0,0,0
        ];
        let s = Silhouette::new(points, "open".to_string());
        assert!(!s.is_closed());
    }

    #[test]
    fn test_zero_area_detection() {
        // Collinear points
        let points = vec![[0.0, 0.0, 0.0], [5.0, 0.0, 0.0], [10.0, 0.0, 0.0]];
        let s = Silhouette::new(points, "flat".to_string());
        assert!(s.get_area() < 1e-6);
    }

    #[test]
    fn test_emerged_and_submerged_equivalence() {
        let rect = create_rectangle(0.0, 100.0, 0.0, 20.0);

        // Fully emerged (waterline at 0.0) vs fully submerged (waterline at 25.0)
        let emerged_area = rect.get_emerged_area(0.0);
        let submerged_area = rect.get_submerged_area(25.0);

        assert!(
            (emerged_area - submerged_area).abs() < 1e-6,
            "Fully emerged area ({}) should equal fully submerged area ({})",
            emerged_area,
            submerged_area
        );

        let emerged_centroid = rect.get_emerged_centroid(0.0);
        let submerged_centroid = rect.get_submerged_centroid(25.0);

        assert!(
            (emerged_centroid[0] - submerged_centroid[0]).abs() < 1e-6
                && (emerged_centroid[1] - submerged_centroid[1]).abs() < 1e-6,
            "Fully emerged centroid ({:?}) should equal fully submerged centroid ({:?})",
            emerged_centroid,
            submerged_centroid
        );
    }

    /// Regression test: submerged centroid of a rectangle fully below the waterline.
    ///
    /// For a 20 × 1.6 rectangle (keel at z=0, top at z=1.6) with waterline at z=1.6
    /// the centroid must be at (10.0, 0.8) — not at the waterline.
    #[test]
    fn test_submerged_centroid_full_rectangle() {
        let rect = create_rectangle(0.0, 20.0, 0.0, 1.6);
        let centroid = rect.get_submerged_centroid(1.6);
        assert!(
            (centroid[0] - 10.0).abs() < 1e-4,
            "Expected x=10.0, got {}",
            centroid[0]
        );
        assert!(
            (centroid[1] - 0.8).abs() < 1e-4,
            "Expected z=0.8, got {}",
            centroid[1]
        );
    }

    /// Regression test: windage-only silhouette (starts at the waterline).
    ///
    /// When the silhouette spans only z ∈ [draft, depth] the submerged area must
    /// be negligibly small (no meaningful part of the silhouette is below the waterline).
    /// Before the fix, clip_below() returned a hair-thin strip whose centroid was
    /// erroneously reported at z ≈ draft instead of 0.
    #[test]
    fn test_submerged_area_windage_only_silhouette() {
        // Silhouette from z=1.6 (waterline) to z=2.0 (deck) — above water only
        let windage = create_rectangle(0.0, 20.0, 1.6, 2.0);
        let submerged_area = windage.get_submerged_area(1.6);
        assert!(
            submerged_area < 1e-3,
            "Windage-only silhouette should have negligible submerged area, got {} m²",
            submerged_area
        );
    }

    /// Regression test: polygon_centroid is correct for a plain CCW rectangle.
    ///
    /// This catches the historical bug where `.abs()` was applied to cx/cz
    /// individually rather than dividing by the signed area.
    #[test]
    fn test_centroid_rectangle_ccw() {
        let rect = create_rectangle(0.0, 20.0, 0.0, 1.6);
        let centroid = rect.get_centroid();
        assert!(
            (centroid[0] - 10.0).abs() < 1e-6,
            "Expected x=10.0, got {}",
            centroid[0]
        );
        assert!(
            (centroid[1] - 0.8).abs() < 1e-6,
            "Expected z=0.8, got {}",
            centroid[1]
        );
    }
}