add covariance interpolation mathspec

docs/anise/explanation/frame-safety.md

@@ -9,6 +9,7 @@ ANISE introduces the concept of **Frame Safety** to eliminate these errors.
 In the SPICE toolkit, frames are often referred to by integer IDs. While ANISE maintains compatibility with NAIF IDs, it treats Frames as rich objects. 
 
 Every Frame in ANISE has an `ephemeris_id` and an `orientation_id`, which store the central object identifier and the reference frame orientation identifier. Each frame may also set properties related to the central object, notably:
+
 - **`mu_km3_s2`**: the gravitational parameter in km^3/s^2
 - **`shape`**: the tri-axial ellipsoid shape of the central object, defined by its semi major and semi minor axes and its polar axis.
 
@@ -30,6 +31,7 @@ Before any computation:
 ANISE represents the relationship between frames as a tree. 
 
 When transforming from Frame A to Frame B:
+
 - **Translation Path**: Finds the common ephemeris ancestor (e.g., the Solar System Barycenter) and sums the vectors along the branches.
 - **Rotation Path**: Finds the common orientation ancestor and composes the Direction Cosine Matrices (DCMs) or Quaternions.
 

docs/anise/reference/almanac.md

@@ -15,12 +15,16 @@ The following is a _small subset_ of all the functions available in the Almanac.
 The `Almanac` supports loading various kernel types. You can load files directly or allow ANISE to "guess" the file type.
 
 ### `load`
+
 Loads any supported ANISE or SPICE file.
+
 - **Rust**: `almanac.load("path/to/file")?`
 - **Python**: `almanac.load("path/to/file")`
 
 ### Specialized Loaders
+
 If you know the file type, you can use specialized loaders for better performance or specific configuration:
+
 - `with_spk(spk)`: Add SPK (ephemeris) data.
 - `with_bpc(bpc)`: Add BPC (high-precision orientation) data.
 - `with_planetary_data(pca)`: Add an ANISE planetary constant data kernel (gravity and tri-axial ellipsoid constants and low-fidelity orientation).
@@ -32,17 +36,23 @@ If you know the file type, you can use specialized loaders for better performanc
 The most common use of the `Almanac` is to transform positions and velocities between frames.
 
 ### `transform_to`
+
 Transforms an `Orbit` (state vector) into a different frame at its own epoch.
+
 - **Parameters**: `orbit`, `target_frame`, `aberration`.
 - **Returns**: A new `Orbit` in the target frame.
 
 ### `translate`
+
 Calculates the relative position and velocity between two frames.
+
 - **Parameters**: `target`, `observer`, `epoch`, `aberration`.
 - **Returns**: A `StateVector` (Position + Velocity).
 
 ### `rotate`
+
 Calculates the rotation (DCM) from one orientation frame to another.
+
 - **Parameters**: `target`, `observer`, `epoch`.
 - **Returns**: A 3x3 Direction Cosine Matrix.
 
@@ -51,10 +61,13 @@ Calculates the rotation (DCM) from one orientation frame to another.
 The `Almanac` also stores physical data for bodies defined in the kernels.
 
 ### `frame_info`
+
 Retrieves a `Frame` object containing metadata for a given NAIF ID.
+
 - **Includes**: Gravitational parameter ($\mu$), body radii, and frame type.
 
 ### `angular_velocity_deg_s`
+
 Returns the angular velocity vector in deg/s of the from_frame wrt to the to_frame.
 
 ## Built-in Constants

docs/anise/reference/analysis-api.md

@@ -11,12 +11,14 @@ Use `report_scalars` to generate time-series data for analysis or plotting.
 Expects a `ReportScalars` object containing a list of expressions and a `StateSpec`.
 
 **Parameters:**
+
 - `report`: A `ReportScalars` instance defined with:
     - `state_spec`: The context (Target, Observer, Frame) to evaluate against.
     - `scalars`: A list of `ScalarExpr` items.
 - `time_series`: A `hifitime::TimeSeries` defining the start, end, and step.
 
 **Returns:**
+
 A dictionary (or map) where keys are Epochs and values are maps of `{ "Alias": Value }`.
 
 ```python
@@ -61,12 +63,15 @@ Use `report_events` to find specific instants (points in time).
 ### `almanac.report_events(state, event, start, end)`
 
 Finds discrete events such as:
+
 - **Orbital Events**: Apoapsis, Periapsis.
 - **Min/Max**: Time of maximum elevation, minimum range.
 - **Crossings**: Time when altitude crosses 100km.
 
 **Returns:**
+
 A list of `EventDetails` objects, each containing:
+
 - `epoch`: The precise time of the event.
 - `value`: The value of the scalar at that time.
 - `orbit`: The full spacecraft state at that time.
@@ -95,12 +100,15 @@ Use `report_event_arcs` to find durations.
 ### `almanac.report_event_arcs(state, event, start, end)`
 
 Finds time intervals where a condition is continuously true.
+
 - **Visibility**: "When can I see the station?"
 - **Battery**: "When am I in sunlight?"
 - **Geometry**: "When is the Earth-Sun-Probe angle < 90?"
 
 **Returns:**
+
 A list of `EventArc` objects, each containing:
+
 - `start`: The `EventDetails` of the start (Rising edge).
 - `end`: The `EventDetails` of the end (Falling edge).
 
@@ -118,11 +126,14 @@ A list of `EventArc` objects, each containing:
 ```
 
 **Returns:**
+
 A list of `EventArc` objects, each containing:
+
 - `start`: The `EventDetails` of the start (Rising edge).
 - `end`: The `EventDetails` of the end (Falling edge).
 
 The [`EventArc`](https://docs.rs/anise/latest/anise/analysis/event/struct.EventArc.html#impl-EventArc) provides the following helpers:
+
 - `duration`: the duration of the event, as a hifitime Duration
 - `start_epoch`, `end_epoch`: the epochs of the start and end of the event
 - `midpoint_epoch`: the half-way point in the event, useful if you need to check for some calculation that happens undoubtedly during the event itself.

docs/anise/reference/azimuth-elevation-range.md

@@ -39,7 +39,9 @@ ANISE follows the algorithm described in **Vallado, Section 4.4.3** for SEZ (Sou
 If an `obstructing_body` is provided, ANISE first performs an ellipsoid intersection test. A line segment is drawn from the transmitter to the receiver. If this segment intersects the body's tri-axial ellipsoid, the `obstructed_by` field is populated, and the calculation reflects the obstruction.
 
 ### 2. SEZ Frame Transformation
+
 The transmitter's state is used to define a local **SEZ frame**:
+
 - **Zenith (Z)**: Normal to the reference ellipsoid at the transmitter's location.
 - **South (S)**: Points toward the South pole of the central body.
 - **East (E)**: Completes the right-handed system ($E = Z \times S$).
@@ -49,6 +51,7 @@ Both the transmitter ($\mathbf{r}_{tx}$) and receiver ($\mathbf{r}_{rx}$) positi
 $$\boldsymbol{\rho}_{sez} = \mathbf{r}_{rx\_sez} - \mathbf{r}_{tx\_sez}$$
 
 ### 4. Angle Computation
+
 - **Range**: $\rho = \|\boldsymbol{\rho}_{sez}\|$
 - **Elevation**: $\phi = \arcsin\left(\frac{\rho_z}{\rho}\right)$
 - **Azimuth**: $A = \operatorname{atan2}(\rho_y, -\rho_x)$

docs/anise/reference/frame.md

@@ -5,6 +5,7 @@ Frames define the coordinate systems used for all calculations in ANISE. Every p
 ## Frame Identification
 
 ANISE is compatible with the standard NAIF ID system. You can refer to frames using:
+
 1. **Integer IDs**: e.g., `399` for Earth.
 2. **Predefined Constants**: e.g., `anise::constants::frames::EARTH_J2000`.
 3. **Names**: In some interfaces (like Python), strings can be used if they have been registered in the `Almanac`.
@@ -19,10 +20,12 @@ let orbit = Orbit::cartesian(x, y, z, vx, vy, vz, epoch, frame);
 ```
 
 If you attempt to perform a math operation between two orbits in different frames, ANISE will either:
+
 - **Automatic Transform**: Some APIs will automatically transform the state to a common frame if the `Almanac` is provided.
 - **Error**: If no `Almanac` is provided to resolve the relationship, ANISE will refuse to perform the operation to prevent physical errors.
 
 ## Metadata
+
 A `Frame` object retrieved via `almanac.frame_info()` contains data retrieved from the loaded datasets (SPK, PCK, BPC):
 
 - **`mu_km3_s2`**: The gravitational parameter ($GM$) of the center body in $km^3/s^2$.
@@ -31,7 +34,9 @@ A `Frame` object retrieved via `almanac.frame_info()` contains data retrieved fr
 - **`orientation_id`**: The NAIF ID used for rotation lookups.
 
 ### Orientation and Phase Angles
+
 For planetocentric (body-fixed) frames, ANISE uses the loaded PCK or BPC data to determine the body's orientation. This is defined by three key phase angles:
+
 1. **Right Ascension (RA)** of the body's north pole.
 2. **Declination (Dec)** of the body's north pole.
 3. **Prime Meridian (W)** location at the given epoch.

docs/anise/reference/instrument.md

@@ -7,19 +7,25 @@ Instrument Kernels allow you to define sensors, cameras, and antennas attached t
 An `Instrument` defines the sensor's geometry relative to the spacecraft bus.
 
 ### Mounting
+
 The mounting is defined by a rigid transformation from the **Spacecraft Body Frame** to the **Instrument Frame**.
+
 - **`q_to_i`**: A quaternion (rotation) describing the orientation.
 - **`offset_i`**: A translation vector (lever arm) from the spacecraft center of mass to the sensor origin.
 
 ### Field of View (FOV)
+
 The FOV defines the sensitive volume of the instrument.
+
 - **Conical**: Circular FOV (e.g., dish antennas, LIDARs). Defined by a `half_angle_deg`.
 - **Rectangular**: Box FOV (e.g., camera sensors). Defined by `x_half_angle_deg` (Width) and `y_half_angle_deg` (Height). Assumes the boresight is +Z.
 
 ## Functionality
 
 ### FOV Margin
+
 ANISE can compute the **FOV Margin**, which is the angle between the target and the nearest FOV boundary.
+
 - **Positive**: Target is inside.
 - **Negative**: Target is outside.
 - **Zero**: Target is exactly on the edge.

docs/anise/reference/location.md

@@ -17,6 +17,7 @@ Real-world ground stations do not have a perfect 0-degree horizon. Mountains or
 ANISE supports **Terrain Masks** to account for this.
 
 A `TerrainMask` is a list of Azimuth/Elevation pairs.
+
 - **Azimuth**: The start of the mask segment.
 - **Elevation**: The minimum elevation required for visibility in that segment.
 

docs/anise/reference/mathspec/interpolation/covariance.md

@@ -0,0 +1,70 @@
+Covariances are matrices of expected values, e.g., the square of standard deviations, and are therefore positive semi-definite (PSD) matrices.
+In the case of the CCSDS OEM files, they represent the uncertainty of a Cartesian position and velocity, in either an inertial frame or an orbit-local frame like the [RIC frame](nyxspace/MathSpec/celestial/coord_systems/#ric-frame).
+
+## Foundation
+
+ANISE 0.9.0 introduced covariance interpolation using a **Log-Euclidean Riemannian Interpolation** method, which might be novel in the astrodynamics community. Unlike standard linear element-wise interpolation, this approach
+respects the geometric manifold of Symmetric Positive Definite (SPD) matrices. This guarantees:
+
+1. **Positive Definiteness:** The interpolated covariance matrix is always mathematically valid (all eigenvalues are strictly positive)
+
+2. **Volume Preservation:** It prevents the artificial "swelling" (determinant increase) of uncertainty that occurs when linearly interpolating between two valid matrices. The interpolation follows the "geodesic" (shortest path) on the curved surface of covariance matrices.
+
+
+## Algorithm
+
+_Note:_ While covariances only need to be PSD, this method only works for positive definite matrices, i.e. none of the eigenvalues of the matrix may be zero due to the use of the natural logarithm.
+
+1. Find the nearest covariance before and after the requested epoch, storing the epochs $t_0$ and $t_1$ of the nearest records
+2. Rotate each of these covariances into the desired frame (e.g. if the user wants the interpolated matrix in the RIC frame, rotate both covariances into the RIC frame before proceeding as they are _stable_ points of the interpolation)
+3. For each matrix:
+    1. Compute the Eigendecomposition of the matrix, building a matrix with the eigenvectors and a diagonal matrix of the associated eigenvalues:
+
+        $$ P = Q \Lambda Q^T $$
+
+    1. Ensure all eigenvalues are positive, else this is not a valid covariance because the matrix is not PSD
+    1. Compute the matrix natural logarithm of diagonal matrix, which is simply the natural log of each diagonal element:
+
+        $$ \Lambda ' = diag(\ln(\lambda_i)) $$
+
+    1. Reconstruct the matrix logarithm by mapping the diagonal log-eigenvalues back using the original eigenvectors:
+
+        $$ \ln(P) = Q \Lambda ' Q^T $$
+
+4. Linearly interpolate in the natural logarithm domain after computing the blending factor for the requested epoch:
+
+    $$ \ln(P_k) = (1 - \alpha) \ln(P) + \alpha \ln(P_1) $$
+
+5. Compute the Eigendecomposition of the interpolated matrix, note that $Q_k$ and $\Lambda_k$ are the new eigenvectors and eigenvalues from the previous step:
+
+    $$ ln(P_k) = Q_k \Lambda_k Q_k^T $$
+
+6. Compute the _exponential_ of the diagonal matrix of eigenvalues to map this covariance back into the original manifold:
+
+    $$ P_k = Q_k~diag(\exp(\lambda_{k_i}))~Q_k^T $$
+
+
+## Scalar example
+
+Let $P(0) = \sigma^2_0 = 1$ and $P(1) = \sigma^2_1 = 100$.
+
+With a basic linear interpolation, we would find the half-way covariance to be:
+
+$$ P(0.5) = \frac{\sigma^2_0 + \sigma^2_1}{2} = 50.5 $$
+
+Which leads to a variance of:
+
+$$ sigma_{0.5} = \sqrt{50.5} \simeq 7.10 $$
+
+At the halfway mark, the standard deviation should not be 70% of the final standard deviation.
+Since uncertainty grows geometrically, the halfway point should be between 1 and 10.
+
+**In other words, the uncertainty swells artificially, creating a larger volume than physically warranted.**
+
+Let's compute the same covariance with the Log-Euclidean interpolation method.
+
+$$ P(0.5) = \exp(0.5 \ln(1) + 0.5 \ln(100)) = 10 $$
+
+Which leads to a standard deviation of $sigma_{0.5} \simeq 3.16$.
+
+When considering that a covariance represents a volume, we note that the Log-Euclidean method interpolates the volume log-linearly, preventing the artificial volume inflation seen in linear interpolation.

docs/anise/reference/mathspec/interpolation/index.md

@@ -1,5 +1,6 @@
 ANISE, like SPICE, heavily relies on interpolation for ephemeris and attitude trajectories. This section details each of the algorithms used in ANISE.
 
+- [Covariance](./covariance.md)
 - [Chebyshev](./chebyshev.md)
 - [Hermite](./hermite.md)
-- [Lagrange](./lagrange.md)
+- [Lagrange](./lagrange.md)

docs/anise/reference/orbit.md

@@ -5,6 +5,7 @@ In ANISE, an `Orbit` represents the full state of an object at a specific time w
 ## The `Orbit` Struct
 
 An `Orbit` consists of:
+
 - **State Vector**: Position and Velocity ($\mathbf{r}, \mathbf{v}$).
 - **Epoch**: The time at which the state is defined (using `hifitime`).
 - **Frame**: The reference frame in which the vectors are expressed.
@@ -33,14 +34,18 @@ Defines a state relative to a body's surface. Useful for ground stations.
 Returns a new `Orbit` expressed in the `target_frame`. This accounts for both translation and rotation of the frames.
 
 ### Orbital Elements
+
 You can extract orbital elements from an `Orbit` object:
+
 - `sma_km()`, `ecc()`, `inc_deg()`, `raan_deg()`, `aop_deg()`, `ta_deg()`
 - **Other Anomalies**: `ma_deg()` (Mean Anomaly), `ea_deg()` (Eccentric Anomaly).
 - **Other Parameters**: `periapsis_km()`, `apoapsis_km()`, `hmag()` (Specific Angular Momentum).
 - **Python**: Accessible as methods (e.g., `orbit.sma_km()`).
 
 
 ### Propagation (Two-Body)
+
 ANISE provides basic two-body propagation for quick estimates:
+
 - `at_epoch(new_epoch)`: Propagates the orbit to a new time using Keplerian motion.
 - **Warning**: For high-fidelity propagation including J2, solar radiation pressure, etc., use a dedicated propagator like those found in the `Nyx` library.