affn 0.7.0

//! # Cartesian Direction (Unit Vectors)
//!
//! This module defines [`Direction<F>`], a **dimensionless unit vector** representing
//! orientation in 3D space.
//!
//! ## Mathematical Model
//!
//! A direction is a unit-length vector representing pure orientation without magnitude:
//!
//! - **Frame-dependent**: Orientation is relative to a reference frame `F`
//! - **Center-independent**: Directions are translation-invariant (free vectors)
//! - **Dimensionless**: No length unit; magnitude is always 1
//!
//! ## Key Properties
//!
//! 1. **Immutable normalization**: All constructors produce unit vectors
//! 2. **No translation**: Directions cannot be "moved" to a different origin
//! 3. **Rotation-only transforms**: Frame transformations are the only valid spatial ops
//!
//! ## Supported Operations
//!
//! | Operation | Result | Meaning |
//! |-----------|--------|---------|
//! | `Direction * Length` | `Vector` | Scale direction to displacement |
//! | `normalize(Vector)` | `Direction` | Extract orientation from displacement |
//! | `Position.direction()` | `Direction` | Unit vector from center to position |
//!
//! ## Forbidden Operations (do not compile)
//!
//! - `Direction + Direction` — Adding unit vectors is geometrically invalid
//! - `Direction + Vector` — Translating a direction is meaningless
//! - Center transformations on Direction — Directions have no center
//!
//! ## Example
//!
//! ```rust
//! use affn::cartesian::Direction;
//! use affn::frames::ReferenceFrame;
//!
//! #[derive(Debug, Copy, Clone)]
//! struct WorldFrame;
//! impl ReferenceFrame for WorldFrame {
//!     fn frame_name() -> &'static str { "WorldFrame" }
//! }
//!
//! // Create a normalized direction
//! let dir = Direction::<WorldFrame>::new(1.0, 2.0, 2.0);
//! assert!((dir.x() - 1.0/3.0).abs() < 1e-12);
//! assert!((dir.y() - 2.0/3.0).abs() < 1e-12);
//! assert!((dir.z() - 2.0/3.0).abs() < 1e-12);
//!
//! // Scale to create a Vector
//! use qtty::units::*; use qtty::{Quantity, M, KM, DEG, RAD, SEC}; use qtty::angular::{Degrees, Radians}; use qtty::length::{Meters, Kilometers};
//! let vec = dir.scale(10.0 * M);
//! ```

use super::vector::Displacement;
use super::xyz::XYZ;
use crate::centers::ReferenceCenter;
use crate::frames::ReferenceFrame;
use qtty::length::LengthUnit;
use qtty::Quantity;

use std::marker::PhantomData;
use std::ops::Mul;

/// A unit vector representing orientation in 3D space.
///
/// Directions are frame-dependent but center-independent (free vectors).
/// The internal storage is a `Vector3<f64>` with magnitude 1 (dimensionless).
///
/// # Type Parameters
/// - `F`: The reference frame (e.g., `ICRS`, `EclipticMeanJ2000`, `Equatorial`)
///
/// # Invariants
///
/// All public constructors ensure the direction is normalized. For unchecked
/// construction, use [`from_xyz_unchecked`](Self::from_xyz_unchecked).
///
/// # Zero-Cost Abstraction
///
/// This type is `#[repr(transparent)]` over `XYZ<f64>`, ensuring no runtime
/// overhead compared to raw `Vector3<f64>`.
///
/// # Renormalization Policy
///
/// The `Mul<Direction>` implementation for [`crate::Rotation3`] uses
/// [`Direction::new_unchecked`] internally and does **not** auto-renormalize
/// the result. A single rotation preserves the unit-norm invariant to within
/// machine epsilon, but long composition chains
/// (`r_n * ... * r_2 * r_1 * dir`) accumulate floating-point error and the
/// magnitude will drift away from `1.0`. For pipelines that apply many
/// rotations in sequence (e.g. integrators, repeated frame transforms),
/// call [`renormalize`](Self::renormalize) or
/// [`renormalized`](Self::renormalized) periodically to restore the
/// invariant.
#[repr(transparent)]
#[derive(Debug, Clone, Copy)]
pub struct Direction<F: ReferenceFrame> {
    pub(in crate::cartesian) xyz: XYZ<f64>,
    _frame: PhantomData<F>,
}

// =============================================================================
// Constructors (Normalizing)
// =============================================================================

impl<F: ReferenceFrame> Direction<F> {
    /// Creates a direction from components, normalizing to unit length.
    ///
    /// # Panics
    /// Panics if the input vector has zero magnitude (cannot normalize a zero vector).
    /// Use [`Direction::try_new`] for fallible construction when the input may be zero.
    ///
    /// # Example
    /// ```rust
    /// use affn::cartesian::Direction;
    /// use affn::frames::ReferenceFrame;
    ///
    /// #[derive(Debug, Copy, Clone)]
    /// struct WorldFrame;
    /// impl ReferenceFrame for WorldFrame {
    ///     fn frame_name() -> &'static str { "WorldFrame" }
    /// }
    ///
    /// let dir = Direction::<WorldFrame>::new(3.0, 4.0, 0.0);
    /// assert!((dir.x() - 0.6).abs() < 1e-12);
    /// assert!((dir.y() - 0.8).abs() < 1e-12);
    /// ```
    #[inline]
    #[must_use]
    pub fn new(x: f64, y: f64, z: f64) -> Self {
        Self::try_new(x, y, z).expect("Cannot create Direction from zero vector")
    }

    /// Attempts to create a direction, returning `None` if the input is zero.
    #[inline]
    #[must_use]
    pub fn try_new(x: f64, y: f64, z: f64) -> Option<Self> {
        XYZ::new(x, y, z)
            .try_normalize()
            .map(Self::from_xyz_unchecked)
    }

    /// Creates a direction from components (alias for `new`).
    ///
    /// Provided for API symmetry with earlier versions.
    #[inline]
    pub fn normalize(x: f64, y: f64, z: f64) -> Self {
        Self::new(x, y, z)
    }

    /// Creates a direction from a `[f64; 3]` array, normalizing to unit length.
    #[inline]
    pub fn from_array(v: [f64; 3]) -> Self {
        Self::new(v[0], v[1], v[2])
    }
}

// =============================================================================
// Unchecked Constructors (for internal use)
// =============================================================================

impl<F: ReferenceFrame> Direction<F> {
    /// Creates a direction from pre-normalized XYZ storage.
    ///
    /// # Safety
    /// The caller must ensure `xyz` is already a unit vector (magnitude ≈ 1).
    /// No normalization is performed.
    #[inline]
    pub(crate) fn from_xyz_unchecked(xyz: XYZ<f64>) -> Self {
        Self {
            xyz,
            _frame: PhantomData,
        }
    }

    /// Creates a direction from raw components without normalization.
    ///
    /// # Safety
    /// The caller must ensure the components form a unit vector.
    #[inline]
    pub fn new_unchecked(x: f64, y: f64, z: f64) -> Self {
        Self::from_xyz_unchecked(XYZ::new(x, y, z))
    }
}

// =============================================================================
// Renormalization
// =============================================================================

impl<F: ReferenceFrame> Direction<F> {
    /// Restores the unit-norm invariant in place.
    ///
    /// Divides each component by the current [`magnitude`](XYZ::magnitude).
    /// If the magnitude is non-finite (`NaN` or infinite) or smaller than
    /// `f64::EPSILON`, the value is left unchanged rather than panicking
    /// or producing `NaN`s.
    ///
    /// Use this after long chains of `Rotation3 * Direction` operations to
    /// counteract accumulated floating-point drift.
    #[inline]
    pub fn renormalize(&mut self) {
        let mag = self.xyz.magnitude();
        if mag.is_finite() && mag > f64::EPSILON {
            self.xyz = self.xyz.scale(1.0 / mag);
        }
    }

    /// By-value variant of [`renormalize`](Self::renormalize).
    ///
    /// Returns a renormalized copy of `self`. If the current magnitude is
    /// non-finite or below `f64::EPSILON`, the original value is returned
    /// unchanged.
    #[inline]
    #[must_use]
    pub fn renormalized(mut self) -> Self {
        self.renormalize();
        self
    }
}

// =============================================================================
// Component Access
// =============================================================================

impl<F: ReferenceFrame> Direction<F> {
    /// Returns the x-component.
    #[inline]
    pub fn x(&self) -> f64 {
        self.xyz.x()
    }

    /// Returns the y-component.
    #[inline]
    pub fn y(&self) -> f64 {
        self.xyz.y()
    }

    /// Returns the z-component.
    #[inline]
    pub fn z(&self) -> f64 {
        self.xyz.z()
    }

    /// Returns the components as a `[f64; 3]` array.
    #[inline]
    pub fn as_array(&self) -> [f64; 3] {
        *self.xyz.as_array()
    }

    /// Reinterprets this direction as belonging to a different reference frame.
    ///
    /// This is a **zero-cost** operation: the unit-vector components are
    /// preserved unchanged; only the compile-time frame tag is replaced.
    ///
    /// Use after applying a mathematical rotation whose result still carries
    /// the original frame tag.
    #[inline]
    pub fn reinterpret_frame<F2: ReferenceFrame>(self) -> Direction<F2> {
        Direction::new_unchecked(self.x(), self.y(), self.z())
    }
}

// =============================================================================
// Scaling: Direction * Length -> Vector
// =============================================================================

impl<F: ReferenceFrame> Direction<F> {
    /// Scales the direction by a length to produce a displacement vector.
    ///
    /// # Example
    /// ```rust
    /// use affn::cartesian::Direction;
    /// use affn::frames::ReferenceFrame;
    /// use qtty::units::*; use qtty::{Quantity, M, KM, DEG, RAD, SEC}; use qtty::angular::{Degrees, Radians}; use qtty::length::{Meters, Kilometers};
    ///
    /// #[derive(Debug, Copy, Clone)]
    /// struct WorldFrame;
    /// impl ReferenceFrame for WorldFrame {
    ///     fn frame_name() -> &'static str { "WorldFrame" }
    /// }
    ///
    /// let dir = Direction::<WorldFrame>::new(1.0, 0.0, 0.0);
    /// let vec = dir.scale(5.0 * M);
    /// assert!((vec.x().value() - 5.0).abs() < 1e-12);
    /// ```
    #[inline]
    pub fn scale<U: LengthUnit>(&self, magnitude: Quantity<U>) -> Displacement<F, U> {
        Displacement::new(
            magnitude * self.x(),
            magnitude * self.y(),
            magnitude * self.z(),
        )
    }

    /// Creates a position at the given distance from the origin in this direction.
    ///
    /// For centers with `Params = ()`, this is a convenience method.
    #[inline]
    pub fn position<C, U>(&self, magnitude: Quantity<U>) -> super::Position<C, F, U>
    where
        C: ReferenceCenter<Params = ()>,
        U: LengthUnit,
    {
        super::Position::new(
            magnitude * self.x(),
            magnitude * self.y(),
            magnitude * self.z(),
        )
    }

    /// Creates a position with explicit center parameters.
    #[inline]
    pub fn position_with_params<C, U>(
        &self,
        center_params: C::Params,
        magnitude: Quantity<U>,
    ) -> super::Position<C, F, U>
    where
        C: ReferenceCenter,
        U: LengthUnit,
    {
        super::Position::new_with_params(
            center_params,
            magnitude * self.x(),
            magnitude * self.y(),
            magnitude * self.z(),
        )
    }
}

// =============================================================================
// Operator: Direction * Quantity<U> -> Vector
// =============================================================================

impl<F: ReferenceFrame, U: LengthUnit> Mul<Quantity<U>> for Direction<F> {
    type Output = Displacement<F, U>;

    #[inline]
    fn mul(self, magnitude: Quantity<U>) -> Self::Output {
        self.scale(magnitude)
    }
}

forward_ref_binop! { impl[F: ReferenceFrame, U: LengthUnit] Mul, mul for Direction<F>, Quantity<U> }

impl<F: ReferenceFrame, U: LengthUnit> Mul<Direction<F>> for Quantity<U> {
    type Output = Displacement<F, U>;

    #[inline]
    fn mul(self, dir: Direction<F>) -> Self::Output {
        dir.scale(self)
    }
}

forward_ref_binop! { impl[F: ReferenceFrame, U: LengthUnit] Mul, mul for Quantity<U>, Direction<F> }

// =============================================================================
// Geometric Operations
// =============================================================================

impl<F: ReferenceFrame> Direction<F> {
    /// Computes the dot product with another direction.
    ///
    /// Returns cosine of the angle between directions (range: -1 to 1).
    #[inline]
    pub fn dot(&self, other: &Self) -> f64 {
        self.xyz.dot(&other.xyz)
    }

    /// Computes the cross product with another direction.
    ///
    /// The result is normalized if non-zero (perpendicular directions).
    #[inline]
    pub fn cross(&self, other: &Self) -> Option<Self> {
        self.xyz
            .cross(&other.xyz)
            .try_normalize()
            .map(Self::from_xyz_unchecked)
    }

    /// Negates the direction (points in opposite direction).
    #[inline]
    pub fn negate(&self) -> Self {
        Self::from_xyz_unchecked(self.xyz.neg())
    }

    /// Returns the angle between this direction and another, in radians.
    #[inline]
    pub fn angle_to(&self, other: &Self) -> f64 {
        // Clamp to handle floating-point errors at extremes
        self.dot(other).clamp(-1.0, 1.0).acos()
    }
}

// =============================================================================
// Spherical Conversion
// =============================================================================

impl<F: ReferenceFrame> Direction<F> {
    /// Converts this Cartesian direction to spherical coordinates.
    ///
    /// Returns a spherical direction with polar (latitude) and azimuth (longitude)
    /// angles in degrees. Angles are canonicalized to:
    /// - polar in `[-90°, +90°]`
    /// - azimuth in `[0°, 360°)`
    pub fn to_spherical(&self) -> crate::spherical::Direction<F> {
        let (polar, azimuth) = crate::spherical::xyz_to_polar_azimuth(self.x(), self.y(), self.z());
        crate::spherical::Direction::<F>::new_unchecked(polar, azimuth)
    }
}

// =============================================================================
// Display
// =============================================================================

impl<F: ReferenceFrame> Direction<F> {
    /// Returns a formatted string representation.
    pub fn display(&self) -> String {
        format!("{self}")
    }
}

impl_quantity_fmt_triplet! {
    impl[F] for Direction<F>
    where { F: ReferenceFrame, },
    fmt_each: {},
    |this, f, FmtOne| {
        write!(f, "Frame: {}, X: ", F::frame_name())?;
        FmtOne::fmt(&this.x(), f)?;
        write!(f, ", Y: ")?;
        FmtOne::fmt(&this.y(), f)?;
        write!(f, ", Z: ")?;
        FmtOne::fmt(&this.z(), f)
    }
}

// =============================================================================
// Tests
// =============================================================================

#[cfg(test)]
mod tests {
    use super::*;
    // Import the derive
    use crate::DeriveReferenceCenter as ReferenceCenter;
    use crate::DeriveReferenceFrame as ReferenceFrame;
    #[allow(unused_imports)]
    use qtty::angular::{Degrees, Radians};
    #[allow(unused_imports)]
    use qtty::length::{Kilometers, Meters};
    use qtty::units::Meter;
    use qtty::Quantity;
    use qtty::M;
    // Define test-specific frame
    #[derive(Debug, Copy, Clone, ReferenceFrame)]
    struct TestFrame;
    #[derive(Debug, Copy, Clone, ReferenceCenter)]
    struct TestCenter;

    #[derive(Clone, Debug, Default, PartialEq)]
    struct TestParams {
        tag: i32,
    }

    #[derive(Debug, Copy, Clone, ReferenceCenter)]
    #[center(params = TestParams)]
    struct ParamCenter;

    #[test]
    fn test_direction_normalization() {
        let dir = Direction::<TestFrame>::new(3.0, 4.0, 0.0);
        let norm = (dir.x() * dir.x() + dir.y() * dir.y() + dir.z() * dir.z()).sqrt();
        assert!((norm - 1.0).abs() < 1e-12);
        assert!((dir.x() - 0.6).abs() < 1e-12);
        assert!((dir.y() - 0.8).abs() < 1e-12);
    }

    #[test]
    fn test_direction_scale() {
        let dir = Direction::<TestFrame>::new(1.0, 0.0, 0.0);
        let vec = dir.scale(Quantity::<Meter>::new(5.0));
        assert!((vec.x().value() - 5.0).abs() < 1e-12);
        assert!(vec.y().value().abs() < 1e-12);
        assert!(vec.z().value().abs() < 1e-12);
    }

    #[test]
    fn test_direction_dot_product() {
        let a = Direction::<TestFrame>::new(1.0, 0.0, 0.0);
        let b = Direction::<TestFrame>::new(0.0, 1.0, 0.0);
        // Perpendicular directions have dot product 0
        assert!(a.dot(&b).abs() < 1e-12);

        // Same direction has dot product 1
        assert!((a.dot(&a) - 1.0).abs() < 1e-12);

        // Opposite directions have dot product -1
        assert!((a.dot(&a.negate()) + 1.0).abs() < 1e-12);
    }

    #[test]
    fn test_direction_cross_product() {
        let x = Direction::<TestFrame>::new(1.0, 0.0, 0.0);
        let y = Direction::<TestFrame>::new(0.0, 1.0, 0.0);
        let z = x.cross(&y).expect("perpendicular directions");
        assert!(z.x().abs() < 1e-12);
        assert!(z.y().abs() < 1e-12);
        assert!((z.z() - 1.0).abs() < 1e-12);
    }

    #[test]
    fn test_direction_try_new_zero() {
        assert!(Direction::<TestFrame>::try_new(0.0, 0.0, 0.0).is_none());
    }

    #[test]
    fn test_direction_angle() {
        let a = Direction::<TestFrame>::new(1.0, 0.0, 0.0);
        let b = Direction::<TestFrame>::new(0.0, 1.0, 0.0);
        let angle = a.angle_to(&b);
        assert!((angle - std::f64::consts::FRAC_PI_2).abs() < 1e-12);
    }

    #[test]
    fn test_direction_helpers_and_accessors() {
        let dir = Direction::<TestFrame>::normalize(0.0, 3.0, 4.0);
        let arr = dir.as_array();
        assert!((arr[0] - 0.0).abs() < 1e-12);
        assert!((arr[1] - 0.6).abs() < 1e-12);
        assert!((arr[2] - 0.8).abs() < 1e-12);

        let from_arr = Direction::<TestFrame>::from_array([0.0, 3.0, 4.0]);
        assert!((from_arr.y() - 0.6).abs() < 1e-12);
        assert!((from_arr.z() - 0.8).abs() < 1e-12);

        let unchecked = Direction::<TestFrame>::new_unchecked(1.0, 0.0, 0.0);
        assert!((unchecked.x() - 1.0).abs() < 1e-12);
        assert!(unchecked.y().abs() < 1e-12);

        let spherical = unchecked.to_spherical();
        assert!((spherical.polar.value()).abs() < 1e-12);
        assert!((spherical.azimuth.value()).abs() < 1e-12);
    }

    #[test]
    fn test_direction_position_helpers() {
        let dir = Direction::<TestFrame>::new(1.0, 0.0, 0.0);
        let pos = dir.position::<TestCenter, Meter>(2.0 * M);
        assert!((pos.x().value() - 2.0).abs() < 1e-12);
        assert!(pos.y().value().abs() < 1e-12);

        let params = TestParams { tag: 7 };
        let pos_params = dir.position_with_params::<ParamCenter, Meter>(params.clone(), 3.0 * M);
        assert_eq!(pos_params.center_params(), &params);
        assert!((pos_params.x().value() - 3.0).abs() < 1e-12);
    }

    #[test]
    fn test_direction_scaling_operator_left() {
        let dir = Direction::<TestFrame>::new(0.0, 1.0, 0.0);
        let disp: Displacement<TestFrame, Meter> = 4.0 * M * dir;
        assert!((disp.y().value() - 4.0).abs() < 1e-12);
        assert!(disp.x().value().abs() < 1e-12);
    }

    #[test]
    fn test_direction_display() {
        let dir = Direction::<TestFrame>::new(1.0, 0.0, 0.0);
        let text = dir.display();
        assert!(text.contains("Frame: TestFrame"));
        assert!(text.contains("X: 1"));

        let text_prec = format!("{:.3}", dir);
        let expected_x = format!("{:.3}", dir.x());
        assert!(text_prec.contains(&format!("X: {expected_x}")));

        let text_exp = format!("{:.2e}", dir);
        let expected_y = format!("{:.2e}", dir.y());
        assert!(text_exp.contains(&format!("Y: {expected_y}")));
    }

    #[test]
    fn direction_has_xyz_layout() {
        assert_eq!(
            core::mem::size_of::<Direction<TestFrame>>(),
            core::mem::size_of::<XYZ<f64>>()
        );
        assert_eq!(
            core::mem::align_of::<Direction<TestFrame>>(),
            core::mem::align_of::<XYZ<f64>>()
        );
    }

    #[test]
    fn test_direction_renormalize_after_rotation_chain() {
        use crate::Rotation3;
        use qtty::angular::Radians;

        // Tiny xorshift64* PRNG, deterministic seed (no external dep).
        struct Xs(u64);
        impl Xs {
            fn next_u64(&mut self) -> u64 {
                let mut x = self.0;
                x ^= x << 13;
                x ^= x >> 7;
                x ^= x << 17;
                self.0 = x;
                x.wrapping_mul(0x2545F4914F6CDD1D)
            }
            fn next_unit(&mut self) -> f64 {
                // Map to [0, 1) using top 53 bits.
                ((self.next_u64() >> 11) as f64) * (1.0 / ((1u64 << 53) as f64))
            }
            fn next_signed(&mut self) -> f64 {
                self.next_unit() * 2.0 - 1.0
            }
        }

        let mag_of =
            |d: &Direction<TestFrame>| (d.x() * d.x() + d.y() * d.y() + d.z() * d.z()).sqrt();

        let mut rng = Xs(0x00C0_FFEE_DEAD_BEEF_u64);

        // Two parallel runs: one without renormalization, one with periodic.
        let initial = Direction::<TestFrame>::new(1.0, 0.0, 0.0);
        let mut drifting = initial;
        let mut periodic = initial;

        const N: usize = 10_000;
        const PERIOD: usize = 100;

        for i in 1..=N {
            // Random axis (non-zero) and small angle.
            let mut ax = rng.next_signed();
            let mut ay = rng.next_signed();
            let mut az = rng.next_signed();
            let an = (ax * ax + ay * ay + az * az).sqrt();
            if an < 1e-12 {
                ax = 1.0;
                ay = 0.0;
                az = 0.0;
            } else {
                ax /= an;
                ay /= an;
                az /= an;
            }
            // Small angle in (-1e-3, 1e-3) rad to keep many ops geometrically meaningful.
            let angle = Radians::new(rng.next_signed() * 1e-3);
            let rot = Rotation3::from_axis_angle([ax, ay, az], angle);

            drifting = rot * drifting;
            periodic = rot * periodic;

            if i % PERIOD == 0 {
                // Sanity: drifting hasn't blown up yet.
                let drift_mag = mag_of(&drifting);
                assert!(
                    (drift_mag - 1.0).abs() < 1e-3,
                    "drift too large at step {i}: {drift_mag}"
                );

                // Renormalize a *copy* of drifting and confirm it returns to ~1.
                let renormed = drifting.renormalized();
                let renormed_mag = mag_of(&renormed);
                assert!(
                    (renormed_mag - 1.0).abs() < 1e-15,
                    "renormalized magnitude not unit at step {i}: {renormed_mag}"
                );

                // Apply periodic renormalization to the periodic chain in place.
                periodic.renormalize();
                let periodic_mag = mag_of(&periodic);
                assert!(
                    (periodic_mag - 1.0).abs() < 1e-15,
                    "periodic chain not unit immediately after renormalize: {periodic_mag}"
                );
            }
        }

        // Final: the periodic chain must remain very close to unit norm.
        let final_mag = mag_of(&periodic);
        assert!(
            (final_mag - 1.0).abs() < 1e-12,
            "periodic-renormalization chain drifted: {final_mag}"
        );
    }

    #[test]
    fn test_direction_renormalize_leaves_degenerate_unchanged() {
        // A direction whose stored components are NaN should be left alone
        // (no panic, no replacement with NaNs from a divide-by-NaN).
        let mut bad = Direction::<TestFrame>::new_unchecked(f64::NAN, 0.0, 0.0);
        bad.renormalize();
        assert!(bad.x().is_nan());
    }
}