affn 0.12.0

//! # Direction‐type for Spherical Coordinates (unit vectors)
//!
//! A **direction** represents a *unit‐length* pointing vector: spherical
//! coordinates in which the radial component is implicitly fixed to 1.
//!
//! ## Mathematical Foundations
//!
//! Directions are **free vectors**: they live in the vector space, not in affine
//! space. Unlike positions, directions are translation-invariant:
//!
//! - Rotating a direction is valid (frame transformation)
//! - "Translating" a direction to a different origin is **meaningless**
//!
//! Therefore, `Direction<F>` has no center parameter, only a frame `F`.
//!
//! ## Storage
//!
//! A spherical direction stores only the two angles (polar, azimuth). The implicit
//! radius is always 1. This contrasts with [`Position`](super::Position), which
//! stores a radial distance in addition to the angles.
//!
//! ## Angle Convention
//!
//! - **polar (θ)**: Latitude / declination / altitude — range `[-90°, +90°]`
//! - **azimuth (φ)**: Longitude / right ascension / azimuth — normalized to `[0°, 360°)`
//!
//! Constructors canonicalize angles to these ranges.
//!
//! ## Example
//!
//! ```rust
//! use affn::spherical::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" }
//! }
//!
//! // Create a spherical direction
//! let dir = Direction::<WorldFrame>::new_unchecked(45.0 * DEG, 30.0 * DEG);
//! ```

use crate::centers::ReferenceCenter;
use crate::frames::ReferenceFrame;
use qtty::angular::{Degrees, Radians};
use qtty::length::LengthUnit;
use qtty::Quantity;

use core::marker::PhantomData;

// Serde implementations in separate module
#[cfg(feature = "serde")]
#[path = "direction_serde.rs"]
mod direction_serde;

/// A spherical direction (unit vector) in a specific reference frame.
///
/// Directions are frame-dependent but center-independent (free vectors).
/// They cannot undergo center transformations, only frame transformations.
///
/// The direction stores only angles; the implicit radius is always 1.
/// For a spherical coordinate with explicit distance, use [`Position`](super::Position).
///
/// # Type Parameters
/// - `F`: The reference frame (defines axis orientation)
///
/// # Invariants
///
/// - `polar` is in `[-90°, +90°]`
/// - `azimuth` is in `[0°, 360°)`
#[derive(Debug, Clone, Copy)]
pub struct Direction<F: ReferenceFrame> {
    /// Polar angle (θ) - latitude, declination, or altitude, in degrees.
    /// Range: `[-90°, +90°]`
    pub polar: Degrees,
    /// Azimuthal angle (φ) - longitude, right ascension, or azimuth, in degrees.
    /// Range: `[0°, 360°)`
    pub azimuth: Degrees,
    _frame: PhantomData<F>,
}

impl<F: ReferenceFrame> Direction<F> {
    /// Creates a direction from raw angle values without canonicalization.
    ///
    /// The caller must ensure angles are within canonical ranges:
    /// - `polar` in `[-90°, +90°]`
    /// - `azimuth` in `[0°, 360°)`
    ///
    /// For frame-specific constructors that canonicalize automatically, see the
    /// `new` method generated by `#[derive(ReferenceFrame)]` for spherical frames
    /// (e.g., `Direction::<ICRS>::new(ra, dec)`).
    pub const fn new_unchecked(polar: Degrees, azimuth: Degrees) -> Self {
        Self {
            polar,
            azimuth,
            _frame: PhantomData,
        }
    }

    /// Promotes this direction to a full Position with the supplied radial magnitude.
    ///
    /// # Type Parameters
    /// - `C`: The center for the resulting position
    /// - `U`: The length unit for the magnitude
    ///
    /// # Example
    /// ```rust
    /// use affn::spherical::Direction;
    /// use affn::frames::ReferenceFrame;
    /// use affn::centers::ReferenceCenter;
    /// 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" }
    /// }
    ///
    /// #[derive(Debug, Copy, Clone)]
    /// struct WorldOrigin;
    /// impl ReferenceCenter for WorldOrigin {
    ///     type Params = ();
    ///     fn center_name() -> &'static str { "WorldOrigin" }
    /// }
    ///
    /// let dir = Direction::<WorldFrame>::new_unchecked(0.0*DEG, 0.0*DEG);
    /// let pos = dir.position::<WorldOrigin, Meter>(1.0*M);
    /// assert_eq!(pos.distance.value(), 1.0);
    /// ```
    #[must_use]
    pub fn position<C, U>(&self, magnitude: Quantity<U>) -> super::Position<C, F, U>
    where
        C: ReferenceCenter<Params = ()>,
        U: LengthUnit,
    {
        super::Position::new_unchecked(self.polar, self.azimuth, magnitude)
    }

    /// Promotes this direction to a Position with explicit center parameters.
    #[must_use]
    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_unchecked_with_params(
            center_params,
            self.polar,
            self.azimuth,
            magnitude,
        )
    }

    /// Calculates the angular separation between this direction and another
    /// using the Vincenty formula (numerically stable).
    pub fn angular_separation(&self, other: &Self) -> Degrees {
        super::angular_separation_impl(self.polar, self.azimuth, other.polar, other.azimuth)
    }

    /// Returns the oriented great-circle bearing from `self` towards `other`.
    ///
    /// The angle is measured in the tangent plane at `self`, from locally
    /// increasing `polar` (latitude, towards the north pole) towards locally
    /// increasing `azimuth` (longitude, east). Defined results are typed radians
    /// in `[0, 2π)`. In an equatorial frame this is position angle north through
    /// east. Unlike [`Self::angular_separation`], this is an orientation, not an
    /// angular distance, and depends on the frame's polar/azimuth conventions.
    ///
    /// Both directions must use the same frame `F`. A frame does not encode an
    /// observation origin; callers must ensure a common origin whenever their
    /// physical interpretation requires one.
    ///
    /// # Undefined and numerical cases
    ///
    /// Returns `None` for coincident or antipodal directions, a starting point
    /// at either pole (no unique local north), or any non-finite input. As with
    /// [`Self::new_unchecked`], finite angles must satisfy the canonical ranges.
    /// A destination at a pole is valid when the starting point is not a pole.
    ///
    /// For latitudes `a`, `b` and longitude difference `d`, the east and north
    /// tangent components are `cos(b) sin(d)` and
    /// `sin(b-a) + 2 sin(a) cos(b) sin²(d/2)`. Their `atan2(east, north)`
    /// gives the bearing. Near antipodality, first reflect the destination
    /// through the origin and then negate both components. This uses the same
    /// stable formula near coincidence and antipodality. Degree differences
    /// and pole complements preserve small offsets before radian conversion.
    ///
    /// No angular tolerance is used: exact coincidence and antipodality give
    /// an exact zero tangent vector after reduction. If both computed tangent
    /// components vanish (including underflow), returns `None`.
    /// Extremely small offsets may be lost to floating-point rounding or
    /// underflow; bearings near singularities are inherently sensitive to input
    /// error. Normalization maps a rounded `2π` to zero.
    ///
    /// # Example
    ///
    /// ```rust
    /// use affn::{prelude::*, spherical::Direction};
    /// use affn::qtty::{angular::Radians, units::Radian, DEG};
    ///
    /// #[derive(Debug, Copy, Clone, ReferenceFrame)]
    /// struct World;
    ///
    /// let a = Direction::<World>::new_unchecked(0.0 * DEG, 350.0 * DEG);
    /// let b = Direction::<World>::new_unchecked(0.0 * DEG, 10.0 * DEG);
    /// let angle = a.position_angle_to(&b);
    /// let east = (90.0 * DEG).to::<Radian>();
    /// assert!((angle.unwrap() - east).abs() < Radians::new(1e-14));
    /// ```
    #[must_use]
    pub fn position_angle_to(&self, other: &Self) -> Option<Radians> {
        let a = self.polar;
        let mut b = other.polar;
        if [a, b, self.azimuth, other.azimuth]
            .iter()
            .any(|angle| !angle.is_finite())
            || a.abs() == Degrees::QUARTER_TURN
        {
            return None;
        }

        // Reflect B through the origin near antipodality, then negate its
        // tangent components. Subtract 180° before subtracting longitudes to
        // preserve residuals even when the direct difference rounds to 180°.
        let mut delta = other.azimuth - self.azimuth;
        let sign =
            if delta.abs() > Degrees::QUARTER_TURN && delta.abs() <= Degrees::QUARTER_TURN * 3.0 {
                b = -b;
                delta = if delta > Degrees::new(0.0) {
                    (other.azimuth - Degrees::HALF_TURN) - self.azimuth
                } else {
                    other.azimuth - (self.azimuth - Degrees::HALF_TURN)
                };
                -1.0
            } else {
                // Preserve small offsets across the longitude boundary as well.
                if delta > Degrees::HALF_TURN {
                    delta = -((Degrees::FULL_TURN - other.azimuth) + self.azimuth);
                } else if delta < -Degrees::HALF_TURN {
                    delta = (Degrees::FULL_TURN - self.azimuth) + other.azimuth;
                }
                1.0
            };

        // Pole complements retain precision when b-a rounds to ±180° or b
        // approaches ±90°. All differences are formed in the stored degrees.
        let difference = b - a;
        let sin_difference = if difference.abs() > Degrees::QUARTER_TURN {
            difference.signum_const()
                * ((Degrees::QUARTER_TURN - a.abs()) + (Degrees::QUARTER_TURN - b.abs())).sin()
        } else {
            difference.sin()
        };
        let cos_b = (Degrees::QUARTER_TURN - b.abs()).sin();
        let half_sin = (delta * 0.5).sin();
        // Keep angular_separation_impl's direct Vincenty calculation separate:
        // sharing this reduction would add bearing-only branches to distance calls.
        let east = sign * cos_b * delta.sin();
        let north = sign * (sin_difference + 2.0 * a.sin() * cos_b * half_sin * half_sin);
        // Exact coincidence/antipodality now give an exact zero tangent vector;
        // the same check handles numerical underflow without an angular cutoff.
        if east == 0.0 && north == 0.0 {
            return None;
        }
        let angle = Radians::new(qtty::Transcendental::atan2(east, north)).normalize();
        // qtty's remainder can round a tiny negative angle up to a full turn.
        Some(if angle >= Radians::FULL_TURN {
            Radians::new(0.0)
        } else {
            angle
        })
    }

    /// Converts to cartesian Direction.
    pub fn to_cartesian(&self) -> crate::cartesian::Direction<F>
    where
        F: ReferenceFrame,
    {
        use qtty::units::Radian;

        let polar_rad = self.polar.to::<Radian>();
        let azimuth_rad = self.azimuth.to::<Radian>();

        let x = azimuth_rad.cos() * polar_rad.cos();
        let y = azimuth_rad.sin() * polar_rad.cos();
        let z = polar_rad.sin();

        crate::cartesian::Direction::<F>::new(x, y, z)
    }

    /// Constructs a spherical direction from a cartesian Direction.
    ///
    /// The resulting angles are canonicalized:
    /// - `polar` in `[-90°, +90°]`
    /// - `azimuth` in `[0°, 360°)`
    pub fn from_cartesian(cart: &crate::cartesian::Direction<F>) -> Self
    where
        F: ReferenceFrame,
    {
        let (polar, azimuth) = super::xyz_to_polar_azimuth(cart.x(), cart.y(), cart.z());
        Self::new_unchecked(polar, azimuth)
    }
}

impl_quantity_fmt_triplet! {
    impl[F] for Direction<F>
    where { F: ReferenceFrame, },
    fmt_each: {},
    |this, f, FmtOne| {
        let (polar_name, azimuth_name, _) =
            F::spherical_names().unwrap_or(("\u{03b8}", "\u{03c6}", "r"));
        write!(f, "Frame: {}, {}: ", F::frame_name(), polar_name)?;
        FmtOne::fmt(&this.polar, f)?;
        write!(f, ", {}: ", azimuth_name)?;
        FmtOne::fmt(&this.azimuth, f)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    // Import the derives and traits
    use crate::{DeriveReferenceCenter as ReferenceCenter, 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::{DEG, M};
    // Define test-specific frame and center
    #[derive(Debug, Copy, Clone, ReferenceFrame)]
    struct TestFrame;
    #[derive(Debug, Copy, Clone, ReferenceCenter)]
    struct TestCenter;

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

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

    #[test]
    fn creates_valid_spherical_direction() {
        let polar = Degrees::new(45.0);
        let azimuth = Degrees::new(90.0);

        let coord = Direction::<TestFrame>::new_unchecked(polar, azimuth);

        assert_eq!(coord.azimuth.value(), 90.0);
        assert_eq!(coord.polar.value(), 45.0);
    }

    #[test]
    fn displays_coordinate_as_string_correctly() {
        let coord = Direction::<TestFrame>::new_unchecked(Degrees::new(30.0), Degrees::new(60.0));
        let output = coord.to_string();
        assert!(output.contains("θ: 30"), "Missing polar angle");
        assert!(output.contains("φ: 60"), "Missing azimuth");
    }

    #[test]
    fn display_respects_format_specifiers() {
        let coord = Direction::<TestFrame>::new_unchecked(Degrees::new(30.0), Degrees::new(60.0));

        let output_prec = format!("{:.2}", coord);
        let expected_polar_prec = format!("{:.2}", coord.polar);
        assert!(output_prec.contains(&format!("θ: {expected_polar_prec}")));

        let output_exp = format!("{:.3e}", coord);
        let expected_az_exp = format!("{:.3e}", coord.azimuth);
        assert!(output_exp.contains(&format!("φ: {expected_az_exp}")));
    }

    #[test]
    fn maintains_high_precision_on_values() {
        let polar = Degrees::new(45.123_456);
        let azimuth = Degrees::new(90.654_321);

        let coord = Direction::<TestFrame>::new_unchecked(polar, azimuth);

        assert!((coord.polar.value() - 45.123_456).abs() < 1e-6);
        assert!((coord.azimuth.value() - 90.654_321).abs() < 1e-6);
    }

    const EPS: f64 = 1e-6;

    #[test]
    fn position_method_promotes_with_given_radius() {
        let dir = Direction::<TestFrame>::new_unchecked(Degrees::new(-30.0), Degrees::new(120.0));
        let pos = dir.position::<TestCenter, Meter>(Quantity::<Meter>::new(2.0));

        // angles are preserved
        assert!(
            (pos.azimuth.value() - 120.0).abs() < EPS,
            "azimuth mismatch: got {}",
            pos.azimuth.value()
        );
        assert!(
            (pos.polar.value() - (-30.0)).abs() < EPS,
            "polar mismatch: got {}",
            pos.polar.value()
        );

        // distance matches the supplied magnitude
        assert!((pos.distance - 2.0 * M).abs() < EPS * M);
    }

    #[test]
    fn angular_separation_identity() {
        let a = Direction::<TestFrame>::new_unchecked(Degrees::new(45.0), Degrees::new(30.0));
        let sep = a.angular_separation(&a);
        assert!(sep.abs().value() < 1e-10, "expected 0°, got {}", sep);
    }

    fn direction(polar: f64, azimuth: f64) -> Direction<TestFrame> {
        Direction::new_unchecked(polar * DEG, azimuth * DEG)
    }

    fn assert_bearing<F: crate::frames::ReferenceFrame>(
        a: Direction<F>,
        b: Direction<F>,
        expected: Degrees,
        tolerance: Radians,
    ) {
        use qtty::units::Radian;
        let angle = a.position_angle_to(&b).expect("defined bearing");
        assert!(angle >= Radians::new(0.0));
        assert!(angle < Radians::new(core::f64::consts::TAU));
        let error = (angle - expected.to::<Radian>()).abs();
        // Compare on the circle, including the normalization boundary.
        let error = error.value().min(core::f64::consts::TAU - error.value());
        assert!(
            Radians::new(error) <= tolerance,
            "got {angle}, expected {expected}"
        );
    }

    #[test]
    fn position_angle_cardinals_and_destination_poles() {
        let tolerance = Radians::new(8.0 * f64::EPSILON);
        for latitude in [-60.0, 0.0, 60.0] {
            assert_bearing(
                direction(latitude, 40.0),
                direction(80.0, 40.0),
                0.0 * DEG,
                tolerance,
            );
            assert_bearing(
                direction(latitude, 40.0),
                direction(-80.0, 40.0),
                180.0 * DEG,
                tolerance,
            );
            for longitude in [0.0, 120.0, 359.0] {
                assert_bearing(
                    direction(latitude, 40.0),
                    direction(90.0, longitude),
                    0.0 * DEG,
                    tolerance,
                );
                assert_bearing(
                    direction(latitude, 40.0),
                    direction(-90.0, longitude),
                    180.0 * DEG,
                    tolerance,
                );
            }
        }
        assert_bearing(
            direction(0.0, 40.0),
            direction(0.0, 70.0),
            90.0 * DEG,
            tolerance,
        );
        assert_bearing(
            direction(0.0, 40.0),
            direction(0.0, 10.0),
            270.0 * DEG,
            tolerance,
        );
    }

    #[test]
    fn position_angle_normalization_and_wraparound() {
        let tolerance = Radians::new(8.0 * f64::EPSILON);
        for (from, to, expected) in [
            (350.0, 10.0, 90.0),
            (10.0, 350.0, 270.0),
            (359.999999, 0.0, 90.0),
            (0.0, 359.999999, 270.0),
        ] {
            assert_bearing(
                direction(0.0, from),
                direction(0.0, to),
                expected * DEG,
                tolerance,
            );
        }
        let north = direction(45.0, 100.0);
        let just_east = direction(60.0, 100.0 + 1e-10);
        let just_west = direction(60.0, 100.0 - 1e-10);
        let east_angle = north.position_angle_to(&just_east).unwrap();
        let west_angle = north.position_angle_to(&just_west).unwrap();
        assert!(east_angle > Radians::new(0.0) && east_angle < Radians::new(1e-10));
        assert!(west_angle > Radians::new(core::f64::consts::TAU - 1e-10));
        assert!(west_angle < Radians::new(core::f64::consts::TAU));
        // Adding this tiny negative atan2 result to TAU rounds to TAU.
        let rounded = direction(0.0, 1e-16)
            .position_angle_to(&direction(45.0, 0.0))
            .unwrap();
        assert!(rounded >= Radians::new(0.0) && rounded < Radians::new(1e-14));
    }

    #[test]
    fn position_angle_exact_singularities() {
        for latitude in [-89.0, -30.0, 0.0, 30.0, 89.0] {
            for longitude in [0.0, 30.0, 180.0, 330.0] {
                let a = direction(latitude, longitude);
                assert!(a.position_angle_to(&a).is_none());
                let b = direction(-latitude, (longitude + 180.0) % 360.0);
                assert!(a.position_angle_to(&b).is_none());
            }
        }
        for latitude in [-90.0, 90.0] {
            for longitude in [0.0, 90.0, 300.0] {
                let pole = direction(latitude, longitude);
                assert!(pole.position_angle_to(&direction(0.0, 0.0)).is_none());
                assert!(pole.position_angle_to(&pole).is_none());
                assert!(pole
                    .position_angle_to(&direction(-latitude, longitude))
                    .is_none());
            }
        }
    }

    #[test]
    fn position_angle_preserves_adjacent_float_offsets() {
        let tolerance = Radians::new(32.0 * f64::EPSILON);
        let next = |x: f64| f64::from_bits(x.to_bits() + 1);
        // Subtracting these latitudes *after* radian conversion loses precision.
        assert_bearing(
            direction(60.0, 40.0),
            direction(next(60.0), 40.0),
            0.0 * DEG,
            tolerance,
        );
        assert_bearing(
            direction(60.0, 40.0),
            direction(60.0, next(40.0)),
            90.0 * DEG,
            tolerance,
        );
        assert_bearing(
            direction(30.0, 0.0),
            direction(-next(30.0), 180.0),
            180.0 * DEG,
            tolerance,
        );
        assert_bearing(
            direction(0.0, 0.0),
            direction(0.0, next(180.0)),
            270.0 * DEG,
            tolerance,
        );
        // The raw longitude difference rounds to 180, but the inputs are not antipodal.
        assert_bearing(
            direction(0.0, 1e-15),
            direction(0.0, 180.0),
            90.0 * DEG,
            tolerance,
        );
        assert_bearing(
            direction(0.0, 180.0),
            direction(0.0, 1e-15),
            270.0 * DEG,
            tolerance,
        );
        // Both coordinates change by one ULP. The local tangent ratios tend
        // to 1/2 near coincidence, and -4*sqrt(3) : -1 near antipodality.
        assert_bearing(
            direction(60.0, 40.0),
            direction(next(60.0), next(40.0)),
            26.565_051_177_077_99 * DEG,
            tolerance,
        );
        assert_bearing(
            direction(30.0, 40.0),
            direction(-next(30.0), next(220.0)),
            261.786_789_298_261_8 * DEG,
            tolerance,
        );
        let below_pole = f64::from_bits(90.0f64.to_bits() - 1);
        // For two points at the same latitude approaching the north pole,
        // a 10° longitude step has bearing 90° - 10°/2 = 85° in the limit.
        assert_bearing(
            direction(below_pole, 20.0),
            direction(below_pole, 30.0),
            85.0 * DEG,
            tolerance,
        );
        for latitude in [below_pole, -below_pole] {
            assert_bearing(
                direction(latitude, 40.0),
                direction(0.0, 40.0),
                if latitude > 0.0 { 180.0 } else { 0.0 } * DEG,
                tolerance,
            );
            // Opposite meridians, at equal near-polar latitudes: north/south.
            assert_bearing(
                direction(latitude, 40.0),
                direction(latitude, 220.0),
                if latitude > 0.0 { 0.0 } else { 180.0 } * DEG,
                tolerance,
            );
            assert_bearing(
                direction(latitude, 40.0),
                direction(-latitude, 40.0),
                if latitude > 0.0 { 180.0 } else { 0.0 } * DEG,
                tolerance,
            );
        }
    }

    #[test]
    fn position_angle_pole_complements_with_unequal_offsets() {
        let near_pole = f64::from_bits(90.0f64.to_bits() - 1);
        let farther_from_pole = f64::from_bits(90.0f64.to_bits() - 2);
        let tolerance = Radians::new(32.0 * f64::EPSILON);
        // The starting polar offset is one ULP and the destination's is two.
        // At a 90° longitude difference, east:north tends to 2:1. Reflecting
        // latitude gives the southern case; negating the endpoint adds 180°.
        let bearing = 63.434_948_822_922_01;
        for hemisphere in [-1.0, 1.0] {
            let a = direction(hemisphere * near_pole, 0.0);
            let expected = if hemisphere > 0.0 {
                bearing
            } else {
                180.0 - bearing
            };
            assert_bearing(
                a,
                direction(hemisphere * farther_from_pole, 90.0),
                expected * DEG,
                tolerance,
            );
            assert_bearing(
                a,
                direction(-hemisphere * farther_from_pole, 270.0),
                (expected + 180.0) * DEG,
                tolerance,
            );
        }
    }

    #[test]
    fn position_angle_analytic_reference_and_hemisphere_reflection() {
        // At A=(30°,0°), B=(60°,90°), projecting B onto A's tangent basis
        // gives east=1/2 and north=3/4, hence tan(P)=2/3.
        let tolerance = Radians::new(16.0 * f64::EPSILON);
        assert_bearing(
            direction(30.0, 0.0),
            direction(60.0, 90.0),
            33.690_067_525_979_785 * DEG,
            tolerance,
        );
        assert_bearing(
            direction(-30.0, 0.0),
            direction(-60.0, 90.0),
            (180.0 - 33.690_067_525_979_785) * DEG,
            tolerance,
        );
    }

    #[test]
    fn position_angle_matches_constructed_great_circle_tangents() {
        // Construct endpoints in Cartesian space from a known initial tangent
        // and arc length. This does not reuse the spherical bearing formula.
        for latitude in [-75.0, -30.0, 0.0, 45.0, 80.0] {
            for longitude in [0.0, 10.0, 170.0, 350.0] {
                let a = direction(latitude, longitude);
                let (sp, cp) = a.polar.sin_cos();
                let (sl, cl) = a.azimuth.sin_cos();
                let origin = [cp * cl, cp * sl, sp];
                let north = [-sp * cl, -sp * sl, cp];
                let east = [-sl, cl, 0.0];
                for bearing in [0.0, 30.0, 90.0, 145.0, 225.0, 359.0] {
                    let (sb, cb) = (bearing * DEG).sin_cos();
                    for distance in [1.0, 40.0, 120.0, 179.0] {
                        let (sd, cd) = (distance * DEG).sin_cos();
                        let xyz: [f64; 3] = core::array::from_fn(|i| {
                            cd * origin[i] + sd * (cb * north[i] + sb * east[i])
                        });
                        let cart =
                            crate::cartesian::Direction::<TestFrame>::new(xyz[0], xyz[1], xyz[2]);
                        let b = Direction::from_cartesian(&cart);
                        // Endpoint conversion errors are amplified by 1/sin(distance),
                        // at most ~58 here; 4096 eps also covers pole conditioning.
                        assert_bearing(a, b, bearing * DEG, Radians::new(4096.0 * f64::EPSILON));
                    }
                }
            }
        }
    }

    #[test]
    fn position_angle_rejects_non_finite_and_unresolved_inputs() {
        for invalid in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
            for a in [direction(invalid, 0.0), direction(0.0, invalid)] {
                let b = direction(30.0, 20.0);
                assert!(a.position_angle_to(&b).is_none());
                assert!(b.position_angle_to(&a).is_none());
            }
        }
        // The smallest degree offset underflows when converted to radians.
        assert!(direction(0.0, 0.0)
            .position_angle_to(&direction(f64::from_bits(1), 0.0))
            .is_none());
    }

    #[test]
    fn position_angle_is_frame_generic() {
        #[derive(Debug, Copy, Clone, ReferenceFrame)]
        struct OtherFrame;
        assert_bearing(
            Direction::<OtherFrame>::new_unchecked(0.0 * DEG, 0.0 * DEG),
            Direction::<OtherFrame>::new_unchecked(0.0 * DEG, 30.0 * DEG),
            90.0 * DEG,
            Radians::new(8.0 * f64::EPSILON),
        );
    }

    // =============================================================================
    // Canonicalization Tests
    // =============================================================================

    #[test]
    fn canonicalizes_azimuth_to_positive_range() {
        use crate::spherical::canonicalize_azimuth;
        assert!((canonicalize_azimuth(Degrees::new(-90.0)).value() - 270.0).abs() < EPS);
        assert!((canonicalize_azimuth(Degrees::new(450.0)).value() - 90.0).abs() < EPS);
        assert!(canonicalize_azimuth(Degrees::new(-720.0)).value().abs() < EPS);
    }

    #[test]
    fn folds_polar_to_valid_range() {
        use crate::spherical::canonicalize_polar;
        assert!((canonicalize_polar(Degrees::new(100.0)).value() - 80.0).abs() < EPS);
        assert!((canonicalize_polar(Degrees::new(-100.0)).value() - (-80.0)).abs() < EPS);
    }

    #[test]
    fn canonicalizes_angle_pair_across_poles() {
        use crate::spherical::canonicalize_polar_azimuth;

        let (polar, azimuth) = canonicalize_polar_azimuth(100.0 * DEG, 30.0 * DEG);
        assert!((polar.value() - 80.0).abs() < EPS);
        assert!((azimuth.value() - 210.0).abs() < EPS);

        let (polar, azimuth) = canonicalize_polar_azimuth(-100.0 * DEG, 30.0 * DEG);
        assert!((polar.value() + 80.0).abs() < EPS);
        assert!((azimuth.value() - 210.0).abs() < EPS);
    }

    // =============================================================================
    // Roundtrip Tests
    // =============================================================================

    #[test]
    fn roundtrip_spherical_cartesian_direction() {
        let original =
            Direction::<TestFrame>::new_unchecked(Degrees::new(45.0), Degrees::new(30.0));
        let cartesian = original.to_cartesian();
        let recovered = Direction::from_cartesian(&cartesian);

        assert!(
            (recovered.polar.value() - original.polar.value()).abs() < EPS,
            "polar mismatch: {} vs {}",
            recovered.polar.value(),
            original.polar.value()
        );
        assert!(
            (recovered.azimuth.value() - original.azimuth.value()).abs() < EPS,
            "azimuth mismatch: {} vs {}",
            recovered.azimuth.value(),
            original.azimuth.value()
        );
    }

    #[test]
    fn roundtrip_at_poles() {
        // North pole
        let north = Direction::<TestFrame>::new_unchecked(Degrees::new(90.0), Degrees::new(0.0));
        let cart_n = north.to_cartesian();
        let recovered_n = Direction::from_cartesian(&cart_n);
        assert!((recovered_n.polar.value() - 90.0).abs() < EPS);

        // South pole
        let south = Direction::<TestFrame>::new_unchecked(Degrees::new(-90.0), Degrees::new(0.0));
        let cart_s = south.to_cartesian();
        let recovered_s = Direction::from_cartesian(&cart_s);
        assert!((recovered_s.polar.value() - (-90.0)).abs() < EPS);
    }

    #[test]
    fn roundtrip_at_azimuth_boundaries() {
        // At azimuth = 0
        let dir0 = Direction::<TestFrame>::new_unchecked(Degrees::new(30.0), Degrees::new(0.0));
        let cart0 = dir0.to_cartesian();
        let rec0 = Direction::from_cartesian(&cart0);
        assert!((rec0.azimuth.value() - 0.0).abs() < EPS);

        // Near azimuth = 360 (should wrap to near 0)
        let dir360 = Direction::<TestFrame>::new_unchecked(Degrees::new(30.0), Degrees::new(359.9));
        let cart360 = dir360.to_cartesian();
        let rec360 = Direction::from_cartesian(&cart360);
        assert!((rec360.azimuth.value() - 359.9).abs() < EPS);
    }

    #[test]
    fn uses_raw_construction() {
        let dir = Direction::<TestFrame>::new_unchecked(Degrees::new(10.0), Degrees::new(20.0));
        assert!((dir.polar.value() - 10.0).abs() < EPS);
        assert!((dir.azimuth.value() - 20.0).abs() < EPS);
    }

    #[test]
    fn position_with_params_preserves_center() {
        let dir = Direction::<TestFrame>::new_unchecked(Degrees::new(15.0), Degrees::new(25.0));
        let params = TestParams { id: 7 };
        let pos = dir.position_with_params::<ParamCenter, Meter>(params.clone(), 2.5 * M);
        assert_eq!(pos.center_params(), &params);
        assert!((pos.distance - 2.5 * M).abs() < EPS * M);
    }
}