Skip to content

Commit 8e2868f

Browse files
Merge pull request #39 from nyx-space/doc-covariance-interpolation
add covariance interpolation mathspec
2 parents c59b9a8 + c92e119 commit 8e2868f

15 files changed

Lines changed: 2169 additions & 2 deletions

File tree

docs/anise/explanation/frame-safety.md

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@ ANISE introduces the concept of **Frame Safety** to eliminate these errors.
99
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.
1010

1111
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:
12+
1213
- **`mu_km3_s2`**: the gravitational parameter in km^3/s^2
1314
- **`shape`**: the tri-axial ellipsoid shape of the central object, defined by its semi major and semi minor axes and its polar axis.
1415

@@ -30,6 +31,7 @@ Before any computation:
3031
ANISE represents the relationship between frames as a tree.
3132

3233
When transforming from Frame A to Frame B:
34+
3335
- **Translation Path**: Finds the common ephemeris ancestor (e.g., the Solar System Barycenter) and sums the vectors along the branches.
3436
- **Rotation Path**: Finds the common orientation ancestor and composes the Direction Cosine Matrices (DCMs) or Quaternions.
3537

docs/anise/reference/almanac.md

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,12 +15,16 @@ The following is a _small subset_ of all the functions available in the Almanac.
1515
The `Almanac` supports loading various kernel types. You can load files directly or allow ANISE to "guess" the file type.
1616

1717
### `load`
18+
1819
Loads any supported ANISE or SPICE file.
20+
1921
- **Rust**: `almanac.load("path/to/file")?`
2022
- **Python**: `almanac.load("path/to/file")`
2123

2224
### Specialized Loaders
25+
2326
If you know the file type, you can use specialized loaders for better performance or specific configuration:
27+
2428
- `with_spk(spk)`: Add SPK (ephemeris) data.
2529
- `with_bpc(bpc)`: Add BPC (high-precision orientation) data.
2630
- `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
3236
The most common use of the `Almanac` is to transform positions and velocities between frames.
3337

3438
### `transform_to`
39+
3540
Transforms an `Orbit` (state vector) into a different frame at its own epoch.
41+
3642
- **Parameters**: `orbit`, `target_frame`, `aberration`.
3743
- **Returns**: A new `Orbit` in the target frame.
3844

3945
### `translate`
46+
4047
Calculates the relative position and velocity between two frames.
48+
4149
- **Parameters**: `target`, `observer`, `epoch`, `aberration`.
4250
- **Returns**: A `StateVector` (Position + Velocity).
4351

4452
### `rotate`
53+
4554
Calculates the rotation (DCM) from one orientation frame to another.
55+
4656
- **Parameters**: `target`, `observer`, `epoch`.
4757
- **Returns**: A 3x3 Direction Cosine Matrix.
4858

@@ -51,10 +61,13 @@ Calculates the rotation (DCM) from one orientation frame to another.
5161
The `Almanac` also stores physical data for bodies defined in the kernels.
5262

5363
### `frame_info`
64+
5465
Retrieves a `Frame` object containing metadata for a given NAIF ID.
66+
5567
- **Includes**: Gravitational parameter ($\mu$), body radii, and frame type.
5668

5769
### `angular_velocity_deg_s`
70+
5871
Returns the angular velocity vector in deg/s of the from_frame wrt to the to_frame.
5972

6073
## Built-in Constants

docs/anise/reference/analysis-api.md

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,12 +11,14 @@ Use `report_scalars` to generate time-series data for analysis or plotting.
1111
Expects a `ReportScalars` object containing a list of expressions and a `StateSpec`.
1212

1313
**Parameters:**
14+
1415
- `report`: A `ReportScalars` instance defined with:
1516
- `state_spec`: The context (Target, Observer, Frame) to evaluate against.
1617
- `scalars`: A list of `ScalarExpr` items.
1718
- `time_series`: A `hifitime::TimeSeries` defining the start, end, and step.
1819

1920
**Returns:**
21+
2022
A dictionary (or map) where keys are Epochs and values are maps of `{ "Alias": Value }`.
2123

2224
```python
@@ -61,12 +63,15 @@ Use `report_events` to find specific instants (points in time).
6163
### `almanac.report_events(state, event, start, end)`
6264

6365
Finds discrete events such as:
66+
6467
- **Orbital Events**: Apoapsis, Periapsis.
6568
- **Min/Max**: Time of maximum elevation, minimum range.
6669
- **Crossings**: Time when altitude crosses 100km.
6770

6871
**Returns:**
72+
6973
A list of `EventDetails` objects, each containing:
74+
7075
- `epoch`: The precise time of the event.
7176
- `value`: The value of the scalar at that time.
7277
- `orbit`: The full spacecraft state at that time.
@@ -95,12 +100,15 @@ Use `report_event_arcs` to find durations.
95100
### `almanac.report_event_arcs(state, event, start, end)`
96101

97102
Finds time intervals where a condition is continuously true.
103+
98104
- **Visibility**: "When can I see the station?"
99105
- **Battery**: "When am I in sunlight?"
100106
- **Geometry**: "When is the Earth-Sun-Probe angle < 90?"
101107

102108
**Returns:**
109+
103110
A list of `EventArc` objects, each containing:
111+
104112
- `start`: The `EventDetails` of the start (Rising edge).
105113
- `end`: The `EventDetails` of the end (Falling edge).
106114

@@ -118,11 +126,14 @@ A list of `EventArc` objects, each containing:
118126
```
119127

120128
**Returns:**
129+
121130
A list of `EventArc` objects, each containing:
131+
122132
- `start`: The `EventDetails` of the start (Rising edge).
123133
- `end`: The `EventDetails` of the end (Falling edge).
124134

125135
The [`EventArc`](https://docs.rs/anise/latest/anise/analysis/event/struct.EventArc.html#impl-EventArc) provides the following helpers:
136+
126137
- `duration`: the duration of the event, as a hifitime Duration
127138
- `start_epoch`, `end_epoch`: the epochs of the start and end of the event
128139
- `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

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,9 @@ ANISE follows the algorithm described in **Vallado, Section 4.4.3** for SEZ (Sou
3939
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.
4040

4141
### 2. SEZ Frame Transformation
42+
4243
The transmitter's state is used to define a local **SEZ frame**:
44+
4345
- **Zenith (Z)**: Normal to the reference ellipsoid at the transmitter's location.
4446
- **South (S)**: Points toward the South pole of the central body.
4547
- **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
4951
$$\boldsymbol{\rho}_{sez} = \mathbf{r}_{rx\_sez} - \mathbf{r}_{tx\_sez}$$
5052

5153
### 4. Angle Computation
54+
5255
- **Range**: $\rho = \|\boldsymbol{\rho}_{sez}\|$
5356
- **Elevation**: $\phi = \arcsin\left(\frac{\rho_z}{\rho}\right)$
5457
- **Azimuth**: $A = \operatorname{atan2}(\rho_y, -\rho_x)$

docs/anise/reference/frame.md

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ Frames define the coordinate systems used for all calculations in ANISE. Every p
55
## Frame Identification
66

77
ANISE is compatible with the standard NAIF ID system. You can refer to frames using:
8+
89
1. **Integer IDs**: e.g., `399` for Earth.
910
2. **Predefined Constants**: e.g., `anise::constants::frames::EARTH_J2000`.
1011
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);
1920
```
2021

2122
If you attempt to perform a math operation between two orbits in different frames, ANISE will either:
23+
2224
- **Automatic Transform**: Some APIs will automatically transform the state to a common frame if the `Almanac` is provided.
2325
- **Error**: If no `Almanac` is provided to resolve the relationship, ANISE will refuse to perform the operation to prevent physical errors.
2426

2527
## Metadata
28+
2629
A `Frame` object retrieved via `almanac.frame_info()` contains data retrieved from the loaded datasets (SPK, PCK, BPC):
2730

2831
- **`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
3134
- **`orientation_id`**: The NAIF ID used for rotation lookups.
3235

3336
### Orientation and Phase Angles
37+
3438
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:
39+
3540
1. **Right Ascension (RA)** of the body's north pole.
3641
2. **Declination (Dec)** of the body's north pole.
3742
3. **Prime Meridian (W)** location at the given epoch.

docs/anise/reference/instrument.md

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,19 +7,25 @@ Instrument Kernels allow you to define sensors, cameras, and antennas attached t
77
An `Instrument` defines the sensor's geometry relative to the spacecraft bus.
88

99
### Mounting
10+
1011
The mounting is defined by a rigid transformation from the **Spacecraft Body Frame** to the **Instrument Frame**.
12+
1113
- **`q_to_i`**: A quaternion (rotation) describing the orientation.
1214
- **`offset_i`**: A translation vector (lever arm) from the spacecraft center of mass to the sensor origin.
1315

1416
### Field of View (FOV)
17+
1518
The FOV defines the sensitive volume of the instrument.
19+
1620
- **Conical**: Circular FOV (e.g., dish antennas, LIDARs). Defined by a `half_angle_deg`.
1721
- **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.
1822

1923
## Functionality
2024

2125
### FOV Margin
26+
2227
ANISE can compute the **FOV Margin**, which is the angle between the target and the nearest FOV boundary.
28+
2329
- **Positive**: Target is inside.
2430
- **Negative**: Target is outside.
2531
- **Zero**: Target is exactly on the edge.

docs/anise/reference/location.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@ Real-world ground stations do not have a perfect 0-degree horizon. Mountains or
1717
ANISE supports **Terrain Masks** to account for this.
1818

1919
A `TerrainMask` is a list of Azimuth/Elevation pairs.
20+
2021
- **Azimuth**: The start of the mask segment.
2122
- **Elevation**: The minimum elevation required for visibility in that segment.
2223

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
Covariances are matrices of expected values, e.g., the square of standard deviations, and are therefore positive semi-definite (PSD) matrices.
2+
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).
3+
4+
## Foundation
5+
6+
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
7+
respects the geometric manifold of Symmetric Positive Definite (SPD) matrices. This guarantees:
8+
9+
1. **Positive Definiteness:** The interpolated covariance matrix is always mathematically valid (all eigenvalues are strictly positive)
10+
11+
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.
12+
13+
14+
## Algorithm
15+
16+
_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.
17+
18+
1. Find the nearest covariance before and after the requested epoch, storing the epochs $t_0$ and $t_1$ of the nearest records
19+
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)
20+
3. For each matrix:
21+
1. Compute the Eigendecomposition of the matrix, building a matrix with the eigenvectors and a diagonal matrix of the associated eigenvalues:
22+
23+
$$ P = Q \Lambda Q^T $$
24+
25+
1. Ensure all eigenvalues are positive, else this is not a valid covariance because the matrix is not PSD
26+
1. Compute the matrix natural logarithm of diagonal matrix, which is simply the natural log of each diagonal element:
27+
28+
$$ \Lambda ' = diag(\ln(\lambda_i)) $$
29+
30+
1. Reconstruct the matrix logarithm by mapping the diagonal log-eigenvalues back using the original eigenvectors:
31+
32+
$$ \ln(P) = Q \Lambda ' Q^T $$
33+
34+
4. Linearly interpolate in the natural logarithm domain after computing the blending factor for the requested epoch:
35+
36+
$$ \ln(P_k) = (1 - \alpha) \ln(P) + \alpha \ln(P_1) $$
37+
38+
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:
39+
40+
$$ ln(P_k) = Q_k \Lambda_k Q_k^T $$
41+
42+
6. Compute the _exponential_ of the diagonal matrix of eigenvalues to map this covariance back into the original manifold:
43+
44+
$$ P_k = Q_k~diag(\exp(\lambda_{k_i}))~Q_k^T $$
45+
46+
47+
## Scalar example
48+
49+
Let $P(0) = \sigma^2_0 = 1$ and $P(1) = \sigma^2_1 = 100$.
50+
51+
With a basic linear interpolation, we would find the half-way covariance to be:
52+
53+
$$ P(0.5) = \frac{\sigma^2_0 + \sigma^2_1}{2} = 50.5 $$
54+
55+
Which leads to a variance of:
56+
57+
$$ sigma_{0.5} = \sqrt{50.5} \simeq 7.10 $$
58+
59+
At the halfway mark, the standard deviation should not be 70% of the final standard deviation.
60+
Since uncertainty grows geometrically, the halfway point should be between 1 and 10.
61+
62+
**In other words, the uncertainty swells artificially, creating a larger volume than physically warranted.**
63+
64+
Let's compute the same covariance with the Log-Euclidean interpolation method.
65+
66+
$$ P(0.5) = \exp(0.5 \ln(1) + 0.5 \ln(100)) = 10 $$
67+
68+
Which leads to a standard deviation of $sigma_{0.5} \simeq 3.16$.
69+
70+
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.
Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
ANISE, like SPICE, heavily relies on interpolation for ephemeris and attitude trajectories. This section details each of the algorithms used in ANISE.
22

3+
- [Covariance](./covariance.md)
34
- [Chebyshev](./chebyshev.md)
45
- [Hermite](./hermite.md)
5-
- [Lagrange](./lagrange.md)
6+
- [Lagrange](./lagrange.md)

docs/anise/reference/orbit.md

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ In ANISE, an `Orbit` represents the full state of an object at a specific time w
55
## The `Orbit` Struct
66

77
An `Orbit` consists of:
8+
89
- **State Vector**: Position and Velocity ($\mathbf{r}, \mathbf{v}$).
910
- **Epoch**: The time at which the state is defined (using `hifitime`).
1011
- **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.
3334
Returns a new `Orbit` expressed in the `target_frame`. This accounts for both translation and rotation of the frames.
3435

3536
### Orbital Elements
37+
3638
You can extract orbital elements from an `Orbit` object:
39+
3740
- `sma_km()`, `ecc()`, `inc_deg()`, `raan_deg()`, `aop_deg()`, `ta_deg()`
3841
- **Other Anomalies**: `ma_deg()` (Mean Anomaly), `ea_deg()` (Eccentric Anomaly).
3942
- **Other Parameters**: `periapsis_km()`, `apoapsis_km()`, `hmag()` (Specific Angular Momentum).
4043
- **Python**: Accessible as methods (e.g., `orbit.sma_km()`).
4144

4245

4346
### Propagation (Two-Body)
47+
4448
ANISE provides basic two-body propagation for quick estimates:
49+
4550
- `at_epoch(new_epoch)`: Propagates the orbit to a new time using Keplerian motion.
4651
- **Warning**: For high-fidelity propagation including J2, solar radiation pressure, etc., use a dedicated propagator like those found in the `Nyx` library.

0 commit comments

Comments
 (0)