Add frame transformations, geodetic coordinates and topocentric frames

Cartesian::to_frame rotates a state between orientations sharing a center,
velocity included: v' = R v + R_dot r. A test pins it against physics, since
a point fixed to the equator moves at omega*r inertially. ANISE's DCM is
f64, so the matrices are lifted into the propagated scalar; the rotation
depends on epoch and not state, so that is also the right derivative.

Geodetic converts to and from a body-fixed position and gives the local
east-north-up rotation. The height form is the one that survives the poles.

Adds pck11.pca, which reopens body-fixed frames for bodies other than Earth
and lets our body constants be checked against the kernel's. It is fetched
from storage.googleapis.com because public-data.nyxspace.com serves a
certificate that does not cover its own name.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01M9aqCoSmjgoaTbewDHAZhh

data/kernels.tsv

@@ -10,6 +10,12 @@
 # thing a reproducible result can cite, so the dated snapshot covering the same
 # span is pinned instead.
 #
+# pck11.pca comes from ANISE's own bucket via its storage.googleapis.com address
+# rather than public-data.nyxspace.com. That hostname resolves to the same bucket
+# but serves a *.storage.googleapis.com certificate, which does not cover it, so
+# any TLS client refuses it.
+#
 # name	sha256	downloaded	url	purpose
 de440s.bsp	c1c7feeab882263fc493a9d5a5b2ddd71b54826cdf65d8d17a76126b260a49f2	2026-09-30	https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de440s.bsp	Planetary and lunar ephemerides, 1849 to 2150. Sun and Moon positions for third-body and solar radiation pressure.
 earth_000101_261226_260929.bpc	df5510d1995ec56af62c2682affc38781e3aa4fda37722d2fa0a10e531b90b1d	2026-09-30	https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/earth_000101_261226_260929.bpc	Earth orientation, 2000-01-01 to 2026-12-26, reconstructed through 2026-09-29. Needed for inertial to Earth-fixed transformations.
+pck11.pca	4bbe31dd7fa986324343150d0feeaf78a7a102cb9167818b5cd8748db4f3e6a2	2026-09-30	https://storage.googleapis.com/public-data.nyxspace.com/anise/v0.5/pck11.pca	Planetary constants in ANISE's own format: gravitational parameters, shapes and IAU rotation models. Enables body-fixed frames other than Earth's.

src/cartesian.rs

@@ -71,7 +71,9 @@ impl<T: Scalar> Cartesian<T> {
     ///
     /// If the frame is not inertial. Orbital elements describe inertial motion, and
     /// a velocity in a rotating frame is not the inertial velocity — converting one
-    /// anyway would return a plausible-looking wrong answer.
+    /// anyway would return a plausible-looking wrong answer. Call
+    /// [`Cartesian::to_inertial`] first, which needs an epoch and an almanac
+    /// because the rotation is a function of time.
     pub fn keplerian(&self) -> Keplerian<T> {
         Keplerian::from_cartesian(self)
     }

src/frames.rs

@@ -4,9 +4,17 @@
 //! J2000 to ITRF, topocentric, geodetic site coordinates — are a separate piece of
 //! work and will land here on top of ANISE.
 
+use anise::almanac::Almanac;
+use anise::constants::frames::EARTH_ITRF93;
+use anise::constants::orientations::J2000 as J2000_ORIENTATION;
+use hifitime::Epoch;
+use nalgebra::Matrix3;
+use ordinary_diffeq::scalar::Scalar;
 use serde::{Deserialize, Serialize};
+use thiserror::Error;
 
 use crate::body::Body;
+use crate::cartesian::Cartesian;
 
 /// How a frame's axes are oriented.
 ///
@@ -24,9 +32,39 @@ pub enum Orientation {
     BodyFixed,
 }
 
+/// What went wrong transforming between frames.
+#[derive(Debug, Error)]
+pub enum FrameError {
+    /// The two frames are centered on different bodies.
+    #[error(
+        "cannot transform from a frame centered on {from} to one centered on {to}: \
+         only orientation changes are supported, not a change of origin"
+    )]
+    DifferentCenters {
+        /// The origin being transformed from.
+        from: String,
+        /// The origin being transformed to.
+        to: String,
+    },
+
+    /// ANISE could not produce the rotation, usually for want of a kernel
+    /// covering the epoch.
+    #[error("ANISE could not rotate {from} to {to} at {epoch}: {message}")]
+    Rotation {
+        /// The source orientation.
+        from: String,
+        /// The target orientation.
+        to: String,
+        /// The epoch asked for.
+        epoch: Epoch,
+        /// What ANISE said.
+        message: String,
+    },
+}
+
 /// A reference frame: a central body and an orientation.
 ///
-/// Held by every [`Cartesian`](crate::cartesian::Cartesian) and every element set,
+/// Held by every [`Cartesian`] and every element set,
 /// which is what lets conversions read `mu` from the state rather than take it as
 /// an argument.
 #[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
@@ -92,6 +130,110 @@ impl Frame {
             ..*self
         }
     }
+
+    /// The ANISE frame this corresponds to.
+    ///
+    /// Earth's body-fixed frame is **ITRF93**, from the Earth orientation kernel,
+    /// which carries polar motion and the true rotation angle. Every other body's
+    /// is its **IAU** frame, from the planetary constants kernel, which is a
+    /// harmonic model of the pole and prime meridian rather than a measured one —
+    /// good to some tens of meters on the surface, not to centimeters.
+    ///
+    /// Both need their kernel loaded. ANISE reports a missing one when the rotation
+    /// is asked for, which is why this is infallible.
+    pub fn anise(&self) -> anise::prelude::Frame {
+        match self.orientation {
+            Orientation::J2000 => {
+                anise::prelude::Frame::new(self.center.naif_id, J2000_ORIENTATION)
+            }
+            Orientation::BodyFixed if self.center.naif_id == Body::EARTH.naif_id => EARTH_ITRF93,
+            // An IAU body-fixed frame's orientation id is the body's own id.
+            Orientation::BodyFixed => {
+                anise::prelude::Frame::new(self.center.naif_id, self.center.naif_id)
+            }
+        }
+    }
+
+    /// The rotation from this frame's axes to `target`'s, and its time derivative,
+    /// at `epoch`.
+    ///
+    /// Both are real: a frame rotation depends on the epoch, not on the state, so
+    /// there is nothing here for a dual scalar to carry.
+    pub fn rotation_to(
+        &self,
+        target: Frame,
+        epoch: Epoch,
+        almanac: &Almanac,
+    ) -> Result<(Matrix3<f64>, Matrix3<f64>), FrameError> {
+        if self.center.naif_id != target.center.naif_id {
+            return Err(FrameError::DifferentCenters {
+                from: self.center.name(),
+                to: target.center.name(),
+            });
+        }
+
+        if self.orientation == target.orientation {
+            return Ok((Matrix3::identity(), Matrix3::zeros()));
+        }
+
+        let dcm = almanac
+            .rotate(self.anise(), target.anise(), epoch)
+            .map_err(|error| FrameError::Rotation {
+                from: format!("{:?}", self.orientation),
+                to: format!("{:?}", target.orientation),
+                epoch,
+                message: error.to_string(),
+            })?;
+
+        // A rotation without its derivative would drop the frame's own motion from
+        // the velocity, which for Earth-fixed is 465 m/s at the equator.
+        Ok((dcm.rot_mat, dcm.rot_mat_dt.unwrap_or_else(Matrix3::zeros)))
+    }
+}
+
+/// Lifts a real rotation matrix into the propagated scalar.
+///
+/// `nalgebra` is column-major and so is `Matrix3::iter`, so this preserves layout.
+fn lift<T: Scalar>(matrix: &Matrix3<f64>) -> Matrix3<T> {
+    Matrix3::from_iterator(matrix.iter().map(|value| T::lift(*value)))
+}
+
+impl<T: Scalar> Cartesian<T> {
+    /// This state expressed in `target` at `epoch`.
+    ///
+    /// Rotates the position, and the velocity by the same rotation plus the
+    /// frame's own motion: `v' = R v + R_dot r`. Dropping the second term is the
+    /// classic way to get a body-fixed velocity that looks plausible and is
+    /// hundreds of meters per second wrong.
+    ///
+    /// Only orientation changes are supported; the two frames must share a center.
+    pub fn to_frame(
+        &self,
+        target: Frame,
+        epoch: Epoch,
+        almanac: &Almanac,
+    ) -> Result<Self, FrameError> {
+        let (rotation, rotation_rate) = self.frame.rotation_to(target, epoch, almanac)?;
+        let (rotation, rotation_rate) = (lift::<T>(&rotation), lift::<T>(&rotation_rate));
+
+        Ok(Self {
+            position: rotation * self.position,
+            velocity: rotation * self.velocity + rotation_rate * self.position,
+            frame: target,
+        })
+    }
+
+    /// This state in the inertial frame about the same body.
+    ///
+    /// The route from a body-fixed state to orbital elements, which are only
+    /// meaningful for inertial motion.
+    pub fn to_inertial(&self, epoch: Epoch, almanac: &Almanac) -> Result<Self, FrameError> {
+        self.to_frame(
+            self.frame.with_orientation(Orientation::J2000),
+            epoch,
+            almanac,
+        )
+    }
 }
 
 impl Default for Frame {
@@ -103,6 +245,27 @@ impl Default for Frame {
 #[cfg(test)]
 mod tests {
     use super::*;
+    use crate::kernels;
+    use approx::assert_relative_eq;
+    use nalgebra::Vector3;
+    use std::path::{Path, PathBuf};
+
+    /// Earth's rotation rate, rad/s, for checking the frame's own motion.
+    const EARTH_ROTATION_RATE: f64 = 7.292115e-5;
+
+    fn epoch() -> Epoch {
+        Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0)
+    }
+
+    /// The kernels, or `None` when they have not been fetched.
+    fn kernel_directory() -> Option<PathBuf> {
+        let candidate = kernels::directory()
+            .unwrap_or_else(|| Path::new(env!("CARGO_MANIFEST_DIR")).join("data/kernels"));
+        kernels::manifest()
+            .iter()
+            .all(|kernel| candidate.join(kernel.name).is_file())
+            .then_some(candidate)
+    }
 
     #[test]
     fn test_frame_carries_the_central_bodys_mu() {
@@ -125,6 +288,262 @@ mod tests {
         assert_eq!(fixed.center, Body::EARTH);
     }
 
+    #[test]
+    fn test_frames_map_onto_anise() {
+        use anise::constants::frames::IAU_MOON_FRAME;
+
+        // Earth's body-fixed frame is ITRF93, from the orientation kernel.
+        assert_eq!(Frame::EARTH_FIXED.anise(), EARTH_ITRF93);
+        // Every other body's is its IAU frame, from the planetary constants.
+        assert_eq!(Frame::body_fixed(Body::MOON).anise(), IAU_MOON_FRAME);
+        assert_eq!(Frame::MOON_J2000.anise().ephemeris_id, Body::MOON.naif_id);
+    }
+
+    /// A body-fixed rotation about something other than the Earth, which needs the
+    /// planetary constants kernel rather than the Earth orientation one.
+    #[test]
+    fn test_a_non_earth_body_fixed_rotation_works() {
+        let Some(directory) = kernel_directory() else {
+            eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
+            return;
+        };
+        let almanac = kernels::load(&directory).unwrap();
+
+        let (rotation, rate) = Frame::body_fixed(Body::MOON)
+            .rotation_to(Frame::MOON_J2000, epoch(), &almanac)
+            .unwrap();
+        assert_relative_eq!(rotation.determinant(), 1.0, max_relative = 1e-12);
+        // The Moon rotates, so the derivative is not zero.
+        assert!(rate.norm() > 0.0, "rotation rate is zero");
+    }
+
+    /// Our body constants are hardcoded, so they are worth checking against the
+    /// kernel's independently sourced ones.
+    ///
+    /// They are not expected to agree exactly, and the tolerances record by how
+    /// much they do not: our `mu` is EGM96's and the kernel's is DE440's, and our
+    /// shape is WGS84 while the kernel's is the IAU figure. Both differences are
+    /// real, and a run's provenance records which was used.
+    #[test]
+    fn test_body_constants_agree_with_the_kernel() {
+        let Some(directory) = kernel_directory() else {
+            eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
+            return;
+        };
+        let almanac = kernels::load(&directory).unwrap();
+
+        for body in [Body::EARTH, Body::MOON, Body::SUN] {
+            let frame = almanac
+                .frame_info(Frame::inertial(body).anise())
+                .unwrap_or_else(|e| panic!("no constants for {}: {e}", body.name()));
+
+            // EGM96 against DE440 differs in the eighth digit.
+            assert_relative_eq!(
+                body.mu,
+                frame.mu_km3_s2().unwrap() * 1e9,
+                max_relative = 1e-7
+            );
+            // WGS84 against the IAU figure differs by under a meter on Earth.
+            assert_relative_eq!(
+                body.equatorial_radius,
+                frame.mean_equatorial_radius_km().unwrap() * 1e3,
+                max_relative = 1e-6
+            );
+        }
+
+        // Flattening, where Earth's is the one that matters and the Sun is a sphere
+        // in our table but not in the kernel's.
+        assert_relative_eq!(
+            Body::EARTH.flattening,
+            almanac
+                .frame_info(Frame::EARTH_J2000.anise())
+                .unwrap()
+                .flattening()
+                .unwrap(),
+            max_relative = 1e-6
+        );
+    }
+
+    #[test]
+    fn test_a_change_of_center_is_refused() {
+        let Some(directory) = kernel_directory() else {
+            eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
+            return;
+        };
+        let almanac = kernels::load(&directory).unwrap();
+
+        let error = Frame::EARTH_J2000
+            .rotation_to(Frame::MOON_J2000, epoch(), &almanac)
+            .unwrap_err();
+        assert!(error.to_string().contains("change of origin"), "{error}");
+    }
+
+    #[test]
+    fn test_rotation_to_the_same_frame_is_the_identity() {
+        let Some(directory) = kernel_directory() else {
+            eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
+            return;
+        };
+        let almanac = kernels::load(&directory).unwrap();
+
+        let (rotation, rate) = Frame::EARTH_J2000
+            .rotation_to(Frame::EARTH_J2000, epoch(), &almanac)
+            .unwrap();
+        assert_eq!(rotation, Matrix3::identity());
+        assert_eq!(rate, Matrix3::zeros());
+    }
+
+    /// Pins the rotation's *direction* and its derivative against physics rather
+    /// than against itself: a point bolted to the equator is stationary in the
+    /// Earth-fixed frame and moving at `omega * r` in the inertial one.
+    #[test]
+    fn test_a_fixed_point_on_the_equator_moves_at_the_rotation_speed() {
+        let Some(directory) = kernel_directory() else {
+            eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
+            return;
+        };
+        let almanac = kernels::load(&directory).unwrap();
+
+        let radius = Body::EARTH.equatorial_radius;
+        let bolted_down = Cartesian::new(
+            Vector3::new(radius, 0.0, 0.0),
+            Vector3::zeros(),
+            Frame::EARTH_FIXED,
+        );
+
+        let inertial = bolted_down.to_inertial(epoch(), &almanac).unwrap();
+
+        assert_eq!(inertial.frame, Frame::EARTH_J2000);
+        // The rotation preserves length.
+        assert_relative_eq!(inertial.position.norm(), radius, max_relative = 1e-12);
+        // And the frame's own motion appears, to the precision of a nominal
+        // rotation rate against the kernel's actual one.
+        assert_relative_eq!(
+            inertial.velocity.norm(),
+            EARTH_ROTATION_RATE * radius,
+            max_relative = 1e-4
+        );
+        // Moving eastward: perpendicular to the position, and out of the equator's
+        // plane by nothing.
+        assert_relative_eq!(
+            inertial
+                .velocity
+                .normalize()
+                .dot(&inertial.position.normalize()),
+            0.0,
+            epsilon = 1e-9
+        );
+    }
+
+    #[test]
+    fn test_transforming_there_and_back_returns_the_state() {
+        let Some(directory) = kernel_directory() else {
+            eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
+            return;
+        };
+        let almanac = kernels::load(&directory).unwrap();
+
+        let inertial = Cartesian::new(
+            Vector3::new(
+                4.263868426884883e6,
+                5.146189057155391e6,
+                1.1310208421331816e6,
+            ),
+            Vector3::new(-5923.454461876975, 4496.802639690076, 1870.3893008991558),
+            Frame::EARTH_J2000,
+        );
+
+        let fixed = inertial
+            .to_frame(Frame::EARTH_FIXED, epoch(), &almanac)
+            .unwrap();
+        // A LEO velocity in the rotating frame differs by hundreds of m/s, which is
+        // the term a rotation without its derivative would lose.
+        assert!(
+            (fixed.velocity.norm() - inertial.velocity.norm()).abs() > 100.0,
+            "rotating-frame speed {} against inertial {}",
+            fixed.velocity.norm(),
+            inertial.velocity.norm()
+        );
+
+        let back = fixed.to_inertial(epoch(), &almanac).unwrap();
+        assert_relative_eq!(back.position, inertial.position, max_relative = 1e-9);
+        assert_relative_eq!(back.velocity, inertial.velocity, max_relative = 1e-9);
+    }
+
+    /// The path this issue owed: a body-fixed state to orbital elements, which
+    /// `keplerian` alone refuses.
+    #[test]
+    fn test_a_body_fixed_state_reaches_orbital_elements() {
+        let Some(directory) = kernel_directory() else {
+            eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
+            return;
+        };
+        let almanac = kernels::load(&directory).unwrap();
+
+        let inertial = Cartesian::new(
+            Vector3::new(7.0e6, 0.0, 0.0),
+            Vector3::new(0.0, 7546.0, 0.0),
+            Frame::EARTH_J2000,
+        );
+        let fixed = inertial
+            .to_frame(Frame::EARTH_FIXED, epoch(), &almanac)
+            .unwrap();
+
+        let elements = fixed.to_inertial(epoch(), &almanac).unwrap().keplerian();
+        assert_relative_eq!(
+            elements.semi_major_axis,
+            inertial.keplerian().semi_major_axis,
+            max_relative = 1e-9
+        );
+    }
+
+    /// The rotation is state-independent, so it must not contribute to a
+    /// derivative with respect to the state.
+    #[cfg(feature = "autodiff")]
+    #[test]
+    fn test_transforming_preserves_derivatives() {
+        use nalgebra::Vector6;
+        use num_dual::DualSVec64;
+
+        let Some(directory) = kernel_directory() else {
+            eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
+            return;
+        };
+        let almanac = kernels::load(&directory).unwrap();
+
+        // d|r_fixed| / d(state): a rotation preserves length, so the gradient of
+        // the rotated radius is the gradient of the original radius.
+        let state = Vector6::new(7.0e6, 1.0e6, 2.0e6, 0.0, 7546.0, 100.0);
+        let (radius, gradient) = num_dual::gradient(
+            |y: Vector6<DualSVec64<6>>| {
+                Cartesian::new(
+                    y.fixed_rows::<3>(0).into_owned(),
+                    y.fixed_rows::<3>(3).into_owned(),
+                    Frame::EARTH_J2000,
+                )
+                .to_frame(Frame::EARTH_FIXED, epoch(), &almanac)
+                .unwrap()
+                .position
+                .norm()
+            },
+            &state,
+        );
+
+        assert_relative_eq!(
+            radius,
+            state.fixed_rows::<3>(0).norm(),
+            max_relative = 1e-12
+        );
+        let expected = state.fixed_rows::<3>(0).normalize();
+        for i in 0..3 {
+            assert_relative_eq!(gradient[i], expected[i], max_relative = 1e-9);
+        }
+        // Velocity does not enter a position magnitude.
+        for i in 3..6 {
+            assert_relative_eq!(gradient[i], 0.0, epsilon = 1e-12);
+        }
+    }
+
     #[test]
     fn test_frame_round_trips_through_json() {
         let json = serde_json::to_string(&Frame::MOON_J2000).unwrap();

src/geodetic.rs

@@ -0,0 +1,464 @@
+//! Locations on a body's surface, and the local frame at one.
+
+use nalgebra::{Matrix3, Vector3};
+use ordinary_diffeq::scalar::Scalar;
+
+use crate::body::Body;
+use crate::cartesian::Cartesian;
+use crate::frames::{Frame, Orientation};
+
+/// How close to converged the latitude iteration has to get, in radians.
+///
+/// About a nanometer of arc on Earth, which is several orders below anything the
+/// inputs justify.
+const LATITUDE_TOLERANCE: f64 = 1e-14;
+
+/// A location relative to a body's reference ellipsoid.
+///
+/// Latitude is geodetic — measured from the equator to the local vertical, which
+/// on an oblate body is not the direction of the center — and positive north.
+/// Longitude is positive east. Height is along the local vertical from the
+/// ellipsoid, not from the center.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Geodetic<T: Scalar> {
+    /// Geodetic latitude, rad.
+    pub latitude: T,
+    /// Longitude, rad, positive east.
+    pub longitude: T,
+    /// Height above the reference ellipsoid, m.
+    pub height: T,
+    /// The body whose ellipsoid this is measured against.
+    pub body: Body,
+}
+
+impl<T: Scalar> Geodetic<T> {
+    /// A location on `body`.
+    pub fn new(latitude: T, longitude: T, height: T, body: Body) -> Self {
+        Self {
+            latitude,
+            longitude,
+            height,
+            body,
+        }
+    }
+
+    /// The square of the ellipsoid's first eccentricity, `f(2 - f)`.
+    fn eccentricity_squared(&self) -> T {
+        T::lift(self.body.flattening * (2.0 - self.body.flattening))
+    }
+
+    /// Radius of curvature in the prime vertical at this latitude, m.
+    ///
+    /// The distance from the local vertical's intersection with the polar axis out
+    /// to the ellipsoid. Reduces to the radius on a sphere.
+    fn prime_vertical_radius(&self) -> T {
+        let sin_latitude = self.latitude.sin();
+        T::lift(self.body.equatorial_radius)
+            / (T::lift(1.0) - self.eccentricity_squared() * sin_latitude.powi(2)).sqrt()
+    }
+
+    /// This location as a position in the body's body-fixed frame.
+    pub fn to_cartesian(&self) -> Cartesian<T> {
+        let radius = self.prime_vertical_radius();
+        let (sin_latitude, cos_latitude) = (self.latitude.sin(), self.latitude.cos());
+        let (sin_longitude, cos_longitude) = (self.longitude.sin(), self.longitude.cos());
+
+        let equatorial = (radius + self.height) * cos_latitude;
+        Cartesian::new(
+            Vector3::new(
+                equatorial * cos_longitude,
+                equatorial * sin_longitude,
+                (radius * (T::lift(1.0) - self.eccentricity_squared()) + self.height)
+                    * sin_latitude,
+            ),
+            Vector3::zeros(),
+            Frame::body_fixed(self.body),
+        )
+    }
+
+    /// A location from a position in the body's body-fixed frame.
+    ///
+    /// Inverts [`Geodetic::to_cartesian`] by iterating on the latitude, which for
+    /// an Earth-like flattening converges in three or four passes. Under a dual
+    /// scalar the derivative converges with the value.
+    ///
+    /// # Panics
+    ///
+    /// If the state's frame is not body-fixed. Geodetic coordinates rotate with the
+    /// body; reading them off an inertial position would silently give whatever
+    /// longitude the epoch happened to correspond to.
+    pub fn from_cartesian(state: &Cartesian<T>) -> Self {
+        assert!(
+            state.frame.orientation == Orientation::BodyFixed,
+            "geodetic coordinates need a body-fixed frame, but this state is in \
+             {:?} about {}: rotate it first",
+            state.frame.orientation,
+            state.frame.center.name()
+        );
+
+        let body = state.frame.center;
+        let position = state.position;
+        let one = T::lift(1.0);
+        let equatorial = (position[0].powi(2) + position[1].powi(2)).sqrt();
+        let longitude = position[1].atan2(position[0]);
+        let eccentricity_squared = T::lift(body.flattening * (2.0 - body.flattening));
+        let equatorial_radius = T::lift(body.equatorial_radius);
+
+        // The prime vertical radius depends on the latitude and the latitude on it,
+        // so iterate. Start from the geocentric latitude scaled for flattening.
+        let mut latitude = position[2].atan2(equatorial * (one - eccentricity_squared));
+        for _ in 0..10 {
+            let sin_latitude = latitude.sin();
+            let radius =
+                equatorial_radius / (one - eccentricity_squared * sin_latitude.powi(2)).sqrt();
+            let candidate =
+                (position[2] + eccentricity_squared * radius * sin_latitude).atan2(equatorial);
+            let converged = (candidate - latitude).re().abs() < LATITUDE_TOLERANCE;
+            latitude = candidate;
+            if converged {
+                break;
+            }
+        }
+
+        // Height along the local vertical. This form stays well behaved at the poles,
+        // where the more obvious `equatorial / cos(latitude) - radius` does not.
+        let sin_latitude = latitude.sin();
+        let height = equatorial * latitude.cos() + position[2] * sin_latitude
+            - equatorial_radius * (one - eccentricity_squared * sin_latitude.powi(2)).sqrt();
+
+        Self {
+            latitude,
+            longitude,
+            height,
+            body,
+        }
+    }
+
+    /// The rotation from the body-fixed frame to the local east-north-up frame.
+    ///
+    /// Rows are the east, north and up unit vectors in body-fixed components, so
+    /// multiplying a body-fixed vector gives its east, north and up components.
+    /// "Up" is the local vertical, normal to the ellipsoid, which on an oblate body
+    /// does not pass through the center.
+    pub fn east_north_up(&self) -> Matrix3<T> {
+        let (sin_latitude, cos_latitude) = (self.latitude.sin(), self.latitude.cos());
+        let (sin_longitude, cos_longitude) = (self.longitude.sin(), self.longitude.cos());
+        let zero = T::lift(0.0);
+
+        Matrix3::new(
+            -sin_longitude,
+            cos_longitude,
+            zero,
+            -sin_latitude * cos_longitude,
+            -sin_latitude * sin_longitude,
+            cos_latitude,
+            cos_latitude * cos_longitude,
+            cos_latitude * sin_longitude,
+            sin_latitude,
+        )
+    }
+
+    /// Where `target` is, in this location's east-north-up frame.
+    ///
+    /// # Panics
+    ///
+    /// If `target` is not in the same body-fixed frame as this location.
+    pub fn look(&self, target: &Cartesian<T>) -> Vector3<T> {
+        let expected = Frame::body_fixed(self.body);
+        assert!(
+            target.frame == expected,
+            "a look direction needs the target in {:?} about {}, not {:?} about {}",
+            expected.orientation,
+            expected.center.name(),
+            target.frame.orientation,
+            target.frame.center.name()
+        );
+
+        self.east_north_up() * (target.position - self.to_cartesian().position)
+    }
+}
+
+impl<T: Scalar> Cartesian<T> {
+    /// This position as a location on the central body's ellipsoid.
+    ///
+    /// # Panics
+    ///
+    /// If the frame is not body-fixed. Call [`Cartesian::to_frame`] with a
+    /// body-fixed frame first.
+    pub fn geodetic(&self) -> Geodetic<T> {
+        Geodetic::from_cartesian(self)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+    use std::f64::consts::FRAC_PI_2;
+
+    #[test]
+    fn test_equator_on_the_prime_meridian_is_along_x() {
+        let site = Geodetic::new(0.0, 0.0, 0.0, Body::EARTH);
+        let position = site.to_cartesian().position;
+
+        assert_relative_eq!(
+            position[0],
+            Body::EARTH.equatorial_radius,
+            max_relative = 1e-12
+        );
+        assert_relative_eq!(position[1], 0.0, epsilon = 1e-9);
+        assert_relative_eq!(position[2], 0.0, epsilon = 1e-9);
+    }
+
+    /// At the pole the distance from the center is the polar radius, which is where
+    /// flattening shows up most plainly.
+    #[test]
+    fn test_the_pole_is_at_the_polar_radius() {
+        let pole = Geodetic::new(FRAC_PI_2, 0.0, 0.0, Body::EARTH);
+        let position = pole.to_cartesian().position;
+
+        assert_relative_eq!(
+            position[2],
+            Body::EARTH.polar_radius(),
+            max_relative = 1e-12
+        );
+        assert_relative_eq!(
+            position.norm(),
+            Body::EARTH.polar_radius(),
+            max_relative = 1e-12
+        );
+    }
+
+    #[test]
+    fn test_round_trips_over_a_range_of_locations() {
+        let cases = [
+            (0.0_f64, 0.0_f64, 0.0_f64),
+            (45.0, 30.0, 1000.0),
+            (-33.9, 151.2, 58.0),
+            (89.9, -179.9, 0.0),
+            (-89.9, 179.9, 3000.0),
+            (0.0, 180.0, -400.0),
+            (51.5, -0.13, 400_000.0),
+        ];
+
+        for (latitude, longitude, height) in cases {
+            let site = Geodetic::new(
+                latitude.to_radians(),
+                longitude.to_radians(),
+                height,
+                Body::EARTH,
+            );
+            let back = site.to_cartesian().geodetic();
+
+            assert_relative_eq!(
+                back.latitude.to_degrees(),
+                latitude,
+                max_relative = 1e-9,
+                epsilon = 1e-9
+            );
+            assert_relative_eq!(
+                back.longitude.to_degrees(),
+                longitude,
+                max_relative = 1e-9,
+                epsilon = 1e-9
+            );
+            assert_relative_eq!(back.height, height, max_relative = 1e-7, epsilon = 1e-6);
+        }
+    }
+
+    /// Geodetic latitude is not geocentric latitude: on an oblate body the local
+    /// vertical misses the center, by up to about 11.5 arcminutes on Earth.
+    #[test]
+    fn test_geodetic_latitude_differs_from_geocentric() {
+        let site = Geodetic::new(45.0_f64.to_radians(), 0.0, 0.0, Body::EARTH);
+        let position = site.to_cartesian().position;
+
+        let geocentric = position[2].atan2((position[0].powi(2) + position[1].powi(2)).sqrt());
+        let difference = (site.latitude - geocentric).to_degrees() * 60.0;
+        assert!(
+            (11.0..12.0).contains(&difference),
+            "difference {difference} arcminutes, expected about 11.5"
+        );
+    }
+
+    #[test]
+    fn test_east_north_up_is_orthonormal_and_oriented() {
+        let site = Geodetic::new(0.0, 0.0, 0.0, Body::EARTH);
+        let rotation = site.east_north_up();
+
+        assert_relative_eq!(rotation.determinant(), 1.0, max_relative = 1e-12);
+        assert_relative_eq!(
+            rotation * rotation.transpose(),
+            Matrix3::identity(),
+            epsilon = 1e-12
+        );
+
+        // On the equator at longitude zero: up is +x, east is +y, north is +z.
+        assert_relative_eq!(
+            rotation.row(0).transpose(),
+            Vector3::new(0.0, 1.0, 0.0),
+            epsilon = 1e-12
+        );
+        assert_relative_eq!(
+            rotation.row(1).transpose(),
+            Vector3::new(0.0, 0.0, 1.0),
+            epsilon = 1e-12
+        );
+        assert_relative_eq!(
+            rotation.row(2).transpose(),
+            Vector3::new(1.0, 0.0, 0.0),
+            epsilon = 1e-12
+        );
+    }
+
+    /// Up at a site is the local vertical, so a point straight above the site has
+    /// only an up component.
+    #[test]
+    fn test_a_point_overhead_is_purely_up() {
+        let site = Geodetic::new(
+            45.0_f64.to_radians(),
+            30.0_f64.to_radians(),
+            0.0,
+            Body::EARTH,
+        );
+        let overhead = Geodetic::new(site.latitude, site.longitude, 500_000.0, Body::EARTH);
+
+        let look = site.look(&overhead.to_cartesian());
+        assert_relative_eq!(look[0], 0.0, epsilon = 1e-6);
+        assert_relative_eq!(look[1], 0.0, epsilon = 1e-6);
+        assert_relative_eq!(look[2], 500_000.0, max_relative = 1e-9);
+    }
+
+    /// And a point due east along the horizon has no north component.
+    #[test]
+    fn test_east_is_east() {
+        let site = Geodetic::new(0.0, 0.0, 0.0, Body::EARTH);
+        let eastward = Geodetic::new(0.0, 1.0_f64.to_radians(), 0.0, Body::EARTH);
+
+        let look = site.look(&eastward.to_cartesian());
+        assert!(look[0] > 0.0, "east component {}", look[0]);
+        assert_relative_eq!(look[1], 0.0, epsilon = 1e-6);
+        // Below the local horizon, since the surface curves away.
+        assert!(look[2] < 0.0, "up component {}", look[2]);
+    }
+
+    #[test]
+    #[should_panic(expected = "need a body-fixed frame")]
+    fn test_an_inertial_position_is_refused() {
+        let inertial = Cartesian::new(
+            Vector3::new(7.0e6, 0.0, 0.0),
+            Vector3::zeros(),
+            Frame::EARTH_J2000,
+        );
+        let _ = inertial.geodetic();
+    }
+
+    /// A spherical body has no flattening, so geodetic and geocentric latitude
+    /// coincide and the height is just the radius difference.
+    #[test]
+    fn test_a_spherical_body_reduces_to_the_obvious() {
+        let body = Body::spherical(9999, 1.0e14, 1.0e6);
+        let site = Geodetic::new(0.7, 0.3, 1234.0, body);
+        let position = site.to_cartesian().position;
+
+        assert_relative_eq!(position.norm(), 1.0e6 + 1234.0, max_relative = 1e-12);
+        let geocentric = position[2].atan2((position[0].powi(2) + position[1].powi(2)).sqrt());
+        assert_relative_eq!(geocentric, 0.7, max_relative = 1e-12);
+    }
+
+    /// Cross-check against ANISE, whose geodetic conversion is independently
+    /// implemented and validated.
+    #[test]
+    fn test_agrees_with_anise() {
+        use crate::kernels;
+        use anise::constants::frames::EARTH_ITRF93;
+        use anise::prelude::Orbit as AniseOrbit;
+        use hifitime::Epoch;
+        use std::path::{Path, PathBuf};
+
+        let candidate: PathBuf = kernels::directory()
+            .unwrap_or_else(|| Path::new(env!("CARGO_MANIFEST_DIR")).join("data/kernels"));
+        if !kernels::manifest()
+            .iter()
+            .all(|kernel| candidate.join(kernel.name).is_file())
+        {
+            eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
+            return;
+        }
+        let almanac = kernels::load(&candidate).unwrap();
+        let frame = almanac.frame_info(EARTH_ITRF93).unwrap();
+
+        for (latitude, longitude, height) in [(45.0_f64, 30.0_f64, 1000.0), (-33.9, 151.2, 58.0)] {
+            let site = Geodetic::new(
+                latitude.to_radians(),
+                longitude.to_radians(),
+                height,
+                Body::EARTH,
+            );
+            let position = site.to_cartesian().position;
+
+            let theirs = AniseOrbit::new(
+                position[0] * 1e-3,
+                position[1] * 1e-3,
+                position[2] * 1e-3,
+                0.0,
+                0.0,
+                0.0,
+                Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0),
+                frame,
+            );
+
+            // ANISE's ellipsoid comes from the kernel and ours from WGS84, so they
+            // differ by well under a meter rather than exactly.
+            assert_relative_eq!(
+                site.latitude.to_degrees(),
+                theirs.latitude_deg().unwrap(),
+                max_relative = 1e-7
+            );
+            assert_relative_eq!(
+                site.longitude.to_degrees(),
+                theirs.longitude_deg(),
+                max_relative = 1e-9
+            );
+            assert_relative_eq!(
+                site.height,
+                theirs.height_km().unwrap() * 1e3,
+                epsilon = 1.0
+            );
+        }
+    }
+
+    /// The iteration has to differentiate, since a ground station's geometry enters
+    /// a measurement partial.
+    #[cfg(feature = "autodiff")]
+    #[test]
+    fn test_geodetic_conversion_differentiates() {
+        use nalgebra::SVector;
+        use num_dual::DualSVec64;
+
+        let site = Geodetic::new(
+            45.0_f64.to_radians(),
+            30.0_f64.to_radians(),
+            1000.0,
+            Body::EARTH,
+        );
+        let position = site.to_cartesian().position;
+
+        // d(height) / d(position): the local vertical, since height is measured
+        // along it.
+        let (height, gradient): (f64, SVector<f64, 3>) = num_dual::gradient(
+            |p: Vector3<DualSVec64<3>>| {
+                Cartesian::new(p, Vector3::zeros(), Frame::EARTH_FIXED)
+                    .geodetic()
+                    .height
+            },
+            &position,
+        );
+
+        assert_relative_eq!(height, 1000.0, max_relative = 1e-7);
+        let up = site.east_north_up().row(2).transpose();
+        for i in 0..3 {
+            assert_relative_eq!(gradient[i], up[i], max_relative = 1e-6, epsilon = 1e-9);
+        }
+    }
+}

src/kernels.rs

@@ -9,6 +9,11 @@
 //! cites "the latest Earth orientation kernel" cannot be reproduced, so the dated
 //! snapshot covering the same span is pinned instead.
 //!
+//! The planetary constants come from ANISE's bucket by its
+//! `storage.googleapis.com` address. `public-data.nyxspace.com` resolves to the
+//! same bucket but serves a `*.storage.googleapis.com` certificate that does not
+//! cover that name, so a TLS client refuses it.
+//!
 //! A kernel's location is the caller's to decide. [`directory`] reads
 //! `ASTROKIT_KERNELS` for convenience, but nothing here guesses.
 
@@ -141,7 +146,7 @@ mod tests {
     #[test]
     fn test_manifest_parses() {
         let kernels = manifest();
-        assert_eq!(kernels.len(), 2, "{kernels:#?}");
+        assert_eq!(kernels.len(), 3, "{kernels:#?}");
 
         for kernel in &kernels {
             assert!(!kernel.name.is_empty());
@@ -191,6 +196,7 @@ mod tests {
         let almanac = load(&directory).unwrap();
         assert_eq!(almanac.num_loaded_spk(), 1);
         assert_eq!(almanac.num_loaded_bpc(), 1);
+        assert_eq!(almanac.planetary_data.len(), 1);
     }
 
     /// The ephemerides have to actually answer a query, not merely parse.

src/lib.rs

@@ -14,6 +14,7 @@ pub mod cartesian;
 pub mod dynamics;
 pub mod equinoctial;
 pub mod frames;
+pub mod geodetic;
 pub mod kepler;
 pub mod kernels;
 pub mod propagation;
@@ -29,7 +30,8 @@ pub mod prelude {
     pub use crate::cartesian::{Cartesian, FromCartesian};
     pub use crate::dynamics::Dynamics::*;
     pub use crate::equinoctial::{Equinoctial, FromEquinoctial};
-    pub use crate::frames::{Frame, Orientation};
+    pub use crate::frames::{Frame, FrameError, Orientation};
+    pub use crate::geodetic::Geodetic;
     pub use crate::kepler::{Degeneracy, Keplerian, TranslationalLayout};
     pub use crate::propagation::{Accuracy, Propagation};
     pub use crate::state::{