From 432e9cae6532d28a24323d75dc72990cf5dbbf10 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sat, 26 Oct 2024 23:14:55 -0600 Subject: [PATCH 1/6] Update VNC frame --- README.md | 7 ++- .../MathSpec/celestial/coord_systems.md | 45 +++++++++---------- reqs.txt | 6 +++ 3 files changed, 34 insertions(+), 24 deletions(-) create mode 100644 reqs.txt diff --git a/README.md b/README.md index c2ec82e..6d2fa6f 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,12 @@ ## Locally -Run `tasker.sh` to set up the pipenv and run mkdocs, e.g. `tasker.sh serve`. ++ Run `tasker.sh` to set up the pipenv and run mkdocs, e.g. `tasker.sh serve`. ++ Or with [`uv`](https://github.com/astral-sh/uv): + + `uv venv` to create the environment + + `source .venv/bin/activate` to jump into it + + `uv pip install -r reqs.txt` to install the requirements + + `mkdocs serve` to launch the server ### Generating Python documentation diff --git a/docs/nyxspace/MathSpec/celestial/coord_systems.md b/docs/nyxspace/MathSpec/celestial/coord_systems.md index bdb8c1f..0d08491 100644 --- a/docs/nyxspace/MathSpec/celestial/coord_systems.md +++ b/docs/nyxspace/MathSpec/celestial/coord_systems.md @@ -2,46 +2,42 @@ ## Ephemerides -Nyx uses the SPICE developmental ephemerides for all ephemeris handling via ANISE, a modern rewrite of NAIF's SPICE. +Nyx uses the SPICE developmental ephemerides (DE440) for all ephemeris handling via ANISE, a modern rewrite of NAIF's SPICE. Learn all about [ANISE here](../../../anise/). ## Rotations and attitude frames -### Trajectory centered frames -Nyx supports the RIC, VNC and RCN trajectory frames. These frames are right-handed and orthonormal. To retrieve the $3\times 3$ rotation matrix of these frames, call the `dcm_from_traj_frame` function on an `Orbit` structure. -Here is how Nyx computes these frames, where $\Omega$ refers to the [RAAN](orbital_elements.md#right-ascension-of-the-ascending-node-raan), $i$ to the [inclination](orbital_elements.md#inclination-inc), and $u$ to the [argument of latitude](orbital_elements.md#argument-of-latitude-aol). Moreover, $R_1$, $R_2$, $R_3$ respectively correspond to a rotation about the first, second and third axes. +### VNC -#### RIC +The **VNC frame** is a trajectory-centered coordinate system ideal for analyzing spacecraft motion relative to its velocity vector and orbital plane. It is particularly useful for planning maneuvers, such as ensuring a spacecraft burns in the anti-velocity direction during a finite burn. -$$ -[C] = R_3(-\Omega)\times R_1(-i)\times R_3(-u) -$$ +The VNC frame is right-handed and orthonormal, with axes that are mutually perpendicular and unit-length. In some contexts, such as GMAT, it is also referred to as the VNB frame. -#### VNC -We start by computing the unit vectors of the velocity and orbit momentum. +The transformation from the inertial frame to the VNC frame can be represented by a $3 \times 3$ rotation matrix $[C]$. This matrix is constructed using the components of the unit vectors: -$$\mathbf{\hat{v}} = \frac{\mathbf{v}}{v}$$ +1. Velocity, $\mathbf{\hat{v}}$: aligned with the spacecraft's velocity vector in its current frame +2. Normal, $\mathbf{\hat{n}}$: aligned with the orbital momentum, and perpendicular to the orbital plane. +3. Cross-track, $\mathbf{\hat{c}}$: completes the orthonormal basis, perpendicular to both $\mathbf{\hat{v}}$ and $\mathbf{\hat{n}}$. -$$\mathbf{\hat{n}} = \frac{\mathbf{h}}{h}$$ +$$ +[C]= \left[\begin{matrix} +\hat{v}_x & \hat{n}_x & \hat{c}_x \\ +\hat{v}_y & \hat{n}_y & \hat{c}_y \\ +\hat{v}_z & \hat{n}_z & \hat{c}_z \\ +\end{matrix} \right] +$$ -Complete the orthonormal basis: +In practice, you can retrieve this rotation matrix by calling [`dcm_from_vnc_to_inertial`](https://docs.rs/anise/latest/anise/astro/orbit/type.Orbit.html#method.dcm_from_vnc_to_inertial) on an `Orbit` structure, or [`dcm3x3_from_vnc_to_inertial`](https://docs.rs/anise/latest/anise/astro/orbit/type.Orbit.html#method.dcm3x3_from_vnc_to_inertial) for the $3 \times 3$ DCM only. -$$\mathbf{\hat{c}} = \mathbf{\hat{v}} \times \mathbf{\hat{n}}$$ +### RIC $$ -[C]= \left[\begin{matrix} -\hat{v}_x & \hat{n}_x & \hat{c}_x \\ -\hat{v}_y & \hat{n}_y & \hat{c}_y \\ -\hat{v}_z & \hat{n}_z & \hat{c}_z \\ -\end{matrix} -\right] +[C] = R_3(-\Omega)\times R_1(-i)\times R_3(-u) $$ -!!! note - The VNC frame is called VNB in GMAT. -#### RCN +### RCN We start by computing the unit vectors of the radius and orbit momentum. $$\mathbf{\hat{r}} = \frac{\mathbf{r}}{r}$$ @@ -68,4 +64,7 @@ let dcm_vnc2inertial = orbit.dcm_from_traj_frame(Frame::VNC)?; let vector_inertial = dcm_vnc2inertial * vector_vnc; ``` +!!! note + The $6 \times 6$ DCM includes the time derivative of the rotation matrix. ANISE computes this time derivative by a finite difference of the $3 \times 3$ DCMs one millisecond before and after the current epoch. + --8<-- "includes/Abbreviations.md" \ No newline at end of file diff --git a/reqs.txt b/reqs.txt new file mode 100644 index 0000000..8c34d4a --- /dev/null +++ b/reqs.txt @@ -0,0 +1,6 @@ +mkdocs==1.6.1 +mkdocs-git-revision-date-localized-plugin==1.3.0 +mkdocs-glightbox==0.4.0 +mkdocs-jupyter==0.25.1 +mkdocs-material==9.5.42 +mkdocs-redirects==1.2.1 From 8e6b2e3f0de2d5284b3643752b8c7c058d6deebd Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sun, 27 Oct 2024 00:33:26 -0600 Subject: [PATCH 2/6] Add interpolation in ANISE --- .../mathspec/interpolation/chebyshev.md | 75 ++++++++++++++++++ .../mathspec/interpolation/hermite.md | 77 +++++++++++++++++++ .../reference/mathspec/interpolation/index.md | 1 + .../mathspec/interpolation/lagrange.md | 71 +++++++++++++++++ docs/hifitime/index.md | 2 +- .../MathSpec/celestial/coord_systems.md | 68 +++++++++------- mkdocs.yml | 5 ++ 7 files changed, 271 insertions(+), 28 deletions(-) create mode 100644 docs/anise/reference/mathspec/interpolation/chebyshev.md create mode 100644 docs/anise/reference/mathspec/interpolation/hermite.md create mode 100644 docs/anise/reference/mathspec/interpolation/index.md create mode 100644 docs/anise/reference/mathspec/interpolation/lagrange.md diff --git a/docs/anise/reference/mathspec/interpolation/chebyshev.md b/docs/anise/reference/mathspec/interpolation/chebyshev.md new file mode 100644 index 0000000..94bec35 --- /dev/null +++ b/docs/anise/reference/mathspec/interpolation/chebyshev.md @@ -0,0 +1,75 @@ +This implementation efficiently evaluates a Chebyshev polynomial and its derivative using Clenshaw's recurrence formula. The algorithm is optimized for numerical stability and computational efficiency. + +!!! tip + Chebyshev interpolation requires building the Chebyshev polynominal coefficients for each spline of each coordinate of a trajectory. Chebyshev interpolations poorly handle accelerations and as such are not suitable for ephemerides with any kind of eccentricity. They are most commonly used in the interpolation of planetary motion. + +## Foundation + +For Chebyshev polynomials $T_n(x)$ with coefficients $c_n$, the polynomial is: + +$$ P(x) = \frac{c_0}{2} + \sum_{n=1}^N c_n T_n(x) $$ + +where $T_n(x)$ satisfies the recurrence relation: + +$$ T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x) $$ + +## Algorithm + +### Initialization + +- Input: + - Normalized time $x \in [-1,1]$ + - Coefficients $c_n$ + - Spline radius $r_s$ + - Polynomial degree $N$ +- Two three-element working arrays: + - `w`: for function values + - `dw`: for derivative values + +### Clenshaw Recurrence + +For $j = N \to 2$: + +1. Function value recurrence: + + 1. $w_2 = w_1$ + 2. $w_1 = w_0$ + 3. $w_0 = c_j + 2xw_1 - w_2$ + +1. Derivative recurrence: + + 1. $w_2 = w_1$ + 2. $w_1 = w_0$ + 3. $dw_0 = 2w_1 + 2xdw_1 - dw_2$ + +### Final Computation + +1. Function value: + + $$ f(x) = c_0 + xw_0 - w_1 $$ + +2. Derivative (scaled by spline radius): + + $$ f'(x) = \frac{w_0 + xdw_0 - dw_1}{r_s} $$ + +## Error Handling + +The implementation checks for: + +- Non-zero spline radius +- Availability of coefficients at requested indices +- Valid evaluation epoch + +## Complexity + +- Time complexity: $\mathcal{O}(N)$ +- Space complexity: $\mathcal{O}(1)$ (fixed-size working arrays) + +## Key Features + +1. **Numerical Stability**: Uses Clenshaw's algorithm instead of direct polynomial evaluation +2. **Memory Efficiency**: Uses fixed-size arrays of length 3 +3. **Simultaneous Computation**: Evaluates both polynomial and derivative in a single pass +4. **Scale Handling**: Accounts for spline radius in derivative computation + +This implementation is particularly efficient for high-degree polynomials while maintaining numerical stability through the use of Clenshaw recursion. \ No newline at end of file diff --git a/docs/anise/reference/mathspec/interpolation/hermite.md b/docs/anise/reference/mathspec/interpolation/hermite.md new file mode 100644 index 0000000..db1e8e4 --- /dev/null +++ b/docs/anise/reference/mathspec/interpolation/hermite.md @@ -0,0 +1,77 @@ +This implementation computes the Hermite interpolation of a function given a set of points and their derivatives. The algorithm constructs and evaluates both the interpolated function value and its derivative at a specified point. + +!!! tip + Hermite interpolation is particularly useful for trajectory interpolation when provided distinct states with their first order derivative, e.g., position and velocity. + +!!! note + This algorithm differs from the conventional method found in literature as it prevents ringing (aliasing) at specific abscissas. [This test from the original NASA SPICE](https://github.com/nyx-space/anise/blob/0524a6aaa60cca08856260d445b0c53fa6f5c000/anise/src/math/interpolation/hermite.rs#L218) documentation checks that the implementation is free from this ill-behavior. + +## Foundation + +Given $n$ points $(x_i, y_i)$ and their derivatives $y'_i$, the Hermite interpolation constructs a polynomial that matches both the function values and derivatives at each point. The interpolation is built using a divided difference table approach. + +## Algorithm + +[Source code documentation](https://docs.rs/anise/latest/anise/math/interpolation/fn.hermite_eval.html) + +### Initialization + +- Input: Points $(x_i, y_i)$, derivatives $y'_i$, and evaluation point $x$ +- Constraints: $n \leq 32$ points, no duplicate abscissas +- Working array initialized as $2n \times 4n$ matrix + +### First Column Construction + +The initial column alternates between function values and derivatives: + +$$ \text{work}_{2i} = y_i $$ + +$$ \text{work}_{2i+1} = y'_i $$ + +### Second Column Computation + +For each pair of consecutive points $(i-1, i)$: + +1. Compute coefficients: + + 1. $c_1 = x_i - x$ + 1. $c_2 = x - x_{i-1}$ + 1. $\text{denom} = x_i - x_{i-1}$ + +2. Linear interpolation: + + $$ \text{work}_{i-1} = \frac{c_1y_{i-1} + c_2y_i}{\text{denom}} $$ + +### Higher-Order Differences + +For columns $j = 2 \to 2n-1$: + +1. Compute divided differences using: + + $$ \text{work}_{i-1} = \frac{(x_{i+j-1} - x)\text{work}_{i-1} + (x - x_i)\text{work}_i}{x_{i+j-1} - x_i} $$ + +2. Derivative updates: + + $$ \text{work}_{i+2n-1} = \frac{c_1\text{work}_{i+2n-1} + c_2\text{work}_{i+2n} + (\text{work}_i - \text{work}_{i-1})}{\text{denom}} $$ + +## Output + +The algorithm returns a tuple $(f, f')$ where: + +- $f$ is the interpolated function value at $x$ +- $f'$ is the interpolated derivative at $x$ + +## Error Handling + +The implementation includes robust error checking for: + +- Array length consistency +- Maximum sample limit (32 points) +- Division by zero (minimum spacing of $\epsilon \approx 2.22 \times 10^{-16}$) + +## Complexity + +- Time complexity: $\mathcal{O}(n^2)$ +- Space complexity: $\mathcal{O}(n^2)$ + +This algorithm maintains numerical stability through careful handling of divided differences. \ No newline at end of file diff --git a/docs/anise/reference/mathspec/interpolation/index.md b/docs/anise/reference/mathspec/interpolation/index.md new file mode 100644 index 0000000..b077e6f --- /dev/null +++ b/docs/anise/reference/mathspec/interpolation/index.md @@ -0,0 +1 @@ +ANISE, like SPICE, heavily relies on interpolation for ephemeris and attitude trajectories. This section details each of the algorithms used in ANISE. \ No newline at end of file diff --git a/docs/anise/reference/mathspec/interpolation/lagrange.md b/docs/anise/reference/mathspec/interpolation/lagrange.md new file mode 100644 index 0000000..08f813c --- /dev/null +++ b/docs/anise/reference/mathspec/interpolation/lagrange.md @@ -0,0 +1,71 @@ +This implementation computes both the Lagrange interpolation polynomial and its derivative at a specified point using Newton's divided difference method. The algorithm efficiently computes both values simultaneously through a single triangular computation scheme. + +!!! tip + Lagrange interpolation is particularly useful for trajectory interpolation when provided distinct states with their first order derivative, e.g., position and velocity. It is less commonly used in ephemeris interpolation than the [Hermite interpolation](./hermite.md). + +## Foundation + +For $n$ points $(x_i, y_i)$, the Lagrange interpolation polynomial $L(x)$ can be expressed using divided differences. The algorithm computes both: + +$$ L(x) = \sum_{i=0}^{n-1} y_i \prod_{j=0,j\neq i}^{n-1} \frac{x - x_j}{x_i - x_j} $$ + +and its derivative: + +$$ L'(x) = \sum_{i=0}^{n-1} y_i \sum_{k=0,k\neq i}^{n-1} \frac{1}{x_i - x_k} \prod_{j=0,j\neq i,k}^{n-1} \frac{x - x_j}{x_i - x_j} $$ + +## Algorithm + +[Source code documentation](https://docs.rs/anise/latest/anise/math/interpolation/fn.lagrange_eval.html) + +### Initialization + +- Input: Points $(x_i, y_i)$ and evaluation point $x$ +- Two working arrays: + - `work`: stores function values + - `dwork`: stores derivative values + +### Divided Difference Table Construction + +For each level $j = 1 \to n-1$: + +For each index $i = 0 \to n-j-1$: + +1. Compute coefficients: + + $$ x_i, x_{i+j} \text{ (interval endpoints)} $$ + + $$ \text{denom} = x_i - x_{i+j} $$ + +2. Update function value: + + $$ \text{work}_i = \frac{(x - x_{i+j})\text{work}_i + (x_i - x)\text{work}_{i+1}}{\text{denom}} $$ + +3. Update derivative: + + $$ \text{deriv} = \frac{\text{work}_i - \text{work}_{i+1}}{\text{denom}} $$ + + $$ \text{dwork}_i = \frac{(x - x_{i+j})\text{dwork}_i + (x_i - x)\text{dwork}_{i+1}}{\text{denom}} + \text{deriv} $$ + +## Output + +Returns tuple $(f, f')$ where: + +- $f = L(x)$: interpolated function value +- $f' = L'(x)$: interpolated derivative + +## Error Handling + +The implementation includes checks for: + +- Array length consistency +- Non-empty input arrays +- Division by zero (minimum spacing of $\epsilon \approx 2.22 \times 10^{-16}$) + +## Complexity + +The algorithm achieves: + +- Time complexity: $\mathcal{O}(n^2)$ +- Space complexity: $\mathcal{O}(n)$ + +This implementation is particularly efficient as it computes both the interpolated value and its derivative in a single pass through the divided difference table, avoiding redundant calculations. \ No newline at end of file diff --git a/docs/hifitime/index.md b/docs/hifitime/index.md index 39a8aa9..2c99d6e 100644 --- a/docs/hifitime/index.md +++ b/docs/hifitime/index.md @@ -40,7 +40,7 @@ Time scales, or "time systems" as referred to by the European Space Agency, are Humans typically follow Coordinated Universal Time (UTC), which is based on the mean solar time at the prime meridian (0° longitude). Time zones are a longitude and geopolitical construct that apply a fixed offset to UTC such that noon is roughly when the Sun is at its highest point in the sky in a given day. UTC operates on the principle that all seconds are of equal duration, maintaining a constant tick rate over days and years. However, this conflicts with the Earth's rotation, which does not complete in an exact number of seconds and varies unpredictably. -The UT1 time scale measures time relative to the stars, reflecting the Earth's actual rotation. Scientific organizations publish UT1 as part of their daily Earth Orientation Parameters (EOP) updates. Despite these complexities, UTC simplifies matters by assuming a constant 24-hour day, equivalent to 86,400 seconds. +The UT1 time scale measures time relative to the stars, reflecting the Earth's actual rotation. Scientific organizations publish UT1 as part of their daily Earth Orientation Parameters (EOP) updates. Despite these complexities, UTC simplifies matters by assuming a constant 24-hour day, equivalent to 86,400 seconds, unless a leap second is introduced that day. Temps Atomique International (TAI) is a time scale based on the average of atomic clocks worldwide, providing a precise standard for the duration of a second. UTC is an offset from TAI designed to keep it within one second of UT1. This alignment is achieved through leap seconds, which are adjustments announced by the International Earth Rotation and Reference Systems Service (IERS) at least six months in advance. diff --git a/docs/nyxspace/MathSpec/celestial/coord_systems.md b/docs/nyxspace/MathSpec/celestial/coord_systems.md index 0d08491..014a60e 100644 --- a/docs/nyxspace/MathSpec/celestial/coord_systems.md +++ b/docs/nyxspace/MathSpec/celestial/coord_systems.md @@ -1,14 +1,8 @@ # Coordinate systems -## Ephemerides - -Nyx uses the SPICE developmental ephemerides (DE440) for all ephemeris handling via ANISE, a modern rewrite of NAIF's SPICE. - -Learn all about [ANISE here](../../../anise/). - ## Rotations and attitude frames -### VNC +### VNC Frame The **VNC frame** is a trajectory-centered coordinate system ideal for analyzing spacecraft motion relative to its velocity vector and orbital plane. It is particularly useful for planning maneuvers, such as ensuring a spacecraft burns in the anti-velocity direction during a finite burn. @@ -30,41 +24,61 @@ $$ In practice, you can retrieve this rotation matrix by calling [`dcm_from_vnc_to_inertial`](https://docs.rs/anise/latest/anise/astro/orbit/type.Orbit.html#method.dcm_from_vnc_to_inertial) on an `Orbit` structure, or [`dcm3x3_from_vnc_to_inertial`](https://docs.rs/anise/latest/anise/astro/orbit/type.Orbit.html#method.dcm3x3_from_vnc_to_inertial) for the $3 \times 3$ DCM only. -### RIC +Here are the two coordinate frame descriptions following the same style: + +### RIC Frame + +The **RIC frame** is a trajectory-centered coordinate system ideal for analyzing spacecraft motion relative to its position vector and orbital plane. It is particularly useful for relative navigation and rendezvous operations, where spacecraft position relative to a target is critical. + +The RIC frame is right-handed and orthonormal, with axes that are mutually perpendicular and unit-length. + +The transformation from the inertial frame to the RIC frame can be represented by a $3 \times 3$ rotation matrix $[C]$. This matrix is constructed using the components of the unit vectors: + +1. Radial, $\mathbf{\hat{r}}$: aligned with the spacecraft's position vector in its current frame +2. In-track, $\mathbf{\hat{i}}$: perpendicular to both $\mathbf{\hat{r}}$ and $\mathbf{\hat{c}}$, completing the orthonormal basis +3. Cross-track, $\mathbf{\hat{c}}$: aligned with the orbital momentum, perpendicular to the orbital plane $$ -[C] = R_3(-\Omega)\times R_1(-i)\times R_3(-u) +[C]= \left[\begin{matrix} +\hat{r}_x & \hat{i}_x & \hat{c}_x \\ +\hat{r}_y & \hat{i}_y & \hat{c}_y \\ +\hat{r}_z & \hat{i}_z & \hat{c}_z \\ +\end{matrix} \right] $$ +In practice, you can retrieve this rotation matrix by calling [`dcm_from_ric_to_inertial`](https://docs.rs/anise/latest/anise/astro/orbit/type.Orbit.html#method.dcm_from_ric_to_inertial) on an `Orbit` structure, or [`dcm3x3_from_ric_to_inertial`](https://docs.rs/anise/latest/anise/astro/orbit/type.Orbit.html#method.dcm3x3_from_ric_to_inertial) for the $3 \times 3$ DCM only. -### RCN -We start by computing the unit vectors of the radius and orbit momentum. +### RCN Frame -$$\mathbf{\hat{r}} = \frac{\mathbf{r}}{r}$$ +The **RCN frame** is a trajectory-centered coordinate system ideal for analyzing spacecraft motion relative to its position vector and orbital plane. It is particularly useful for closed-loop guidance of finite maneuvers, like in the Q-Law or Ruggiero low-thrust guidance laws. -$$\mathbf{\hat{n}} = \frac{\mathbf{h}}{h}$$ +The RCN frame is right-handed and orthonormal, with axes that are mutually perpendicular and unit-length. -Complete the orthonormal basis: +The transformation from the inertial frame to the RCN frame can be represented by a $3 \times 3$ rotation matrix $[C]$. This matrix is constructed using the components of the unit vectors: -$$\mathbf{\hat{c}} = \mathbf{\hat{r}} \times \mathbf{\hat{n}}$$ +1. Radial, $\mathbf{\hat{r}}$: aligned with the spacecraft's position vector in its current frame +2. Cross-track, $\mathbf{\hat{c}}$: perpendicular to both $\mathbf{\hat{r}}$ and $\mathbf{\hat{n}}$, completing the orthonormal basis +3. Normal, $\mathbf{\hat{n}}$: aligned with the orbital momentum, perpendicular to the orbital plane $$ -[C]= \left[\begin{matrix} -\hat{r}_x & \hat{c}_x & \hat{n}_x \\ -\hat{r}_y & \hat{c}_y & \hat{n}_y \\ -\hat{r}_z & \hat{c}_z & \hat{n}_z \\ -\end{matrix} -\right] +[C]= \left[\begin{matrix} +\hat{r}_x & \hat{c}_x & \hat{n}_x \\ +\hat{r}_y & \hat{c}_y & \hat{n}_y \\ +\hat{r}_z & \hat{c}_z & \hat{n}_z \\ +\end{matrix} \right] $$ -#### Example +In practice, you can retrieve this rotation matrix by calling [`dcm_from_rcn_to_inertial`](https://docs.rs/anise/latest/anise/astro/orbit/type.Orbit.html#method.dcm_from_rcn_to_inertial) on an `Orbit` structure, or [`dcm3x3_from_rcn_to_inertial`](https://docs.rs/anise/latest/anise/astro/orbit/type.Orbit.html#method.dcm3x3_from_rcn_to_inertial) for the $3 \times 3$ DCM only. -``` -let dcm_vnc2inertial = orbit.dcm_from_traj_frame(Frame::VNC)?; -let vector_inertial = dcm_vnc2inertial * vector_vnc; -``` !!! note - The $6 \times 6$ DCM includes the time derivative of the rotation matrix. ANISE computes this time derivative by a finite difference of the $3 \times 3$ DCMs one millisecond before and after the current epoch. + The $6 \times 6$ DCM includes the time derivative of the rotation matrix. ANISE computes this time derivative by a finite difference of the $3 \times 3$ DCMs one millisecond before and after the current epoch. These permutations are handled via [`at_epoch`](https://docs.rs/anise/latest/anise/astro/orbit/type.Orbit.html#method.at_epoch), a two-body approximation of the orbital state using its mean anomaly. + +## Ephemerides + +Nyx uses the NASA/JPL Developmental Ephemerides version 440 (DE440) by default for all ephemeris handling. The implementation is multi-threaded via ANISE, a modern rewrite of NAIF's SPICE. + +Learn all about [ANISE here](../../../anise/). + --8<-- "includes/Abbreviations.md" \ No newline at end of file diff --git a/mkdocs.yml b/mkdocs.yml index 500549f..b14f375 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -105,6 +105,11 @@ nav: - Reference index: anise/reference/index.md - Math Spec: - Orbital elements: anise/reference/mathspec/orbital_elements.md + - Interpolation methods: + - Interpolation index: anise/reference/mathspec/interpolation/index.md + - Chebyshev interpolation: anise/reference/mathspec/interpolation/chebyshev.md + - Hermite interpolation: anise/reference/mathspec/interpolation/hermite.md + - Lagrange interpolation: anise/reference/mathspec/interpolation/lagrange.md - Python API: - ANISE Python: anise/reference/api/python/index.md - Astro: anise/reference/api/python/astro/index.md From d618b070fdffb336fc4c72f6fdae724d5eeed66e Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sun, 27 Oct 2024 01:07:17 -0600 Subject: [PATCH 3/6] Okay, sleep time! --- docs/blog/posts/2024.10.18-hifitime-4.0.0.md | 4 +- .../MathSpec/celestial/coord_systems.md | 2 +- .../MathSpec/celestial/interpolation.md | 143 +++++++++++------- 3 files changed, 92 insertions(+), 57 deletions(-) diff --git a/docs/blog/posts/2024.10.18-hifitime-4.0.0.md b/docs/blog/posts/2024.10.18-hifitime-4.0.0.md index 230cf2d..667df20 100644 --- a/docs/blog/posts/2024.10.18-hifitime-4.0.0.md +++ b/docs/blog/posts/2024.10.18-hifitime-4.0.0.md @@ -15,9 +15,9 @@ links: Months of meticulous development and rigorous testing culminate in version 4.0.0 of Hifitime, a reference time management library for engineering and scientific applications. -Hifitime inherently supports leap seconds, [time scales](/hifitime/#time-scales), Julian dates, and lots more, while guaranteeing exactly one nanosecond accuracy for over 65 thousand years. Hifitime [plans lunar missions](https://fireflyspace.com/missions/blue-ghost-mission-1/) via [ANISE](/anise/), flies onboard various spacecraft, coordinates satellite communications, and is used in exciting research on global navigation satellite systems like GPS and Galileo. +Hifitime inherently supports leap seconds, [time scales](../../hifitime/index.md#time-scales), Julian dates, and lots more, while guaranteeing exactly one nanosecond accuracy for over 65 thousand years. Hifitime [plans lunar missions](https://fireflyspace.com/missions/blue-ghost-mission-1/) via [ANISE](../../anise/index.md), flies onboard various spacecraft, coordinates satellite communications, and is used in exciting research on global navigation satellite systems like GPS and Galileo. -For those eager to dive into the technical details, you can find the [Rust user guide here](/hifitime/rust/) and the [Python user guide here](/hifitime/python/). +For those eager to dive into the technical details, you can find the [Rust user guide here](../../hifitime/rust.md) and the [Python user guide here](../../hifitime/python.md). [Github :material-github:](https://github.com/nyx-space/hifitime){ .md-button } diff --git a/docs/nyxspace/MathSpec/celestial/coord_systems.md b/docs/nyxspace/MathSpec/celestial/coord_systems.md index 014a60e..3a3c6da 100644 --- a/docs/nyxspace/MathSpec/celestial/coord_systems.md +++ b/docs/nyxspace/MathSpec/celestial/coord_systems.md @@ -78,7 +78,7 @@ In practice, you can retrieve this rotation matrix by calling [`dcm_from_rcn_to_ Nyx uses the NASA/JPL Developmental Ephemerides version 440 (DE440) by default for all ephemeris handling. The implementation is multi-threaded via ANISE, a modern rewrite of NAIF's SPICE. -Learn all about [ANISE here](../../../anise/). +Learn all about [ANISE here](../../../anise/index.md). --8<-- "includes/Abbreviations.md" \ No newline at end of file diff --git a/docs/nyxspace/MathSpec/celestial/interpolation.md b/docs/nyxspace/MathSpec/celestial/interpolation.md index ba392f3..b0a4571 100644 --- a/docs/nyxspace/MathSpec/celestial/interpolation.md +++ b/docs/nyxspace/MathSpec/celestial/interpolation.md @@ -1,73 +1,108 @@ -# Interpolated trajectories -This section is focused on the generation of interpolated trajectories by propagators in Nyx. +This section focuses on building interpolated trajectories by propagators in Nyx. -Interpolated trajectories are computed in parallel leading to a negligible impact in computational performance. **All of the event finding requires an interpolated trajectory.** +For details on interpolation methods for trajectories, including [Chebyshev](../../../anise/reference/mathspec/interpolation/chebyshev.md), [Lagrange](../../../anise/reference/mathspec/interpolation/lagrange.md), and [Hermite](../../../anise/reference/mathspec/interpolation/hermite.md) interpolators, refer the [ANISE math spec reference on interpolation](../../../anise/reference/mathspec/interpolation/index.md). -Trajectories store a interpolation splines created from a Hermite interpolation (written from scratch and validated against NumPy). Specifically, the algorithm works by curve fitting the provided states in position, velocity, and fuel mass (when interpolating a spacecraft trajectory instead of just an orbital trajectory). This curve fit will generate a _Hermite series_ which is then converted into a _generic polynomial_. Note that none of these operations require memory allocations (and rely on Rust's [`const generic`](https://rust-lang.github.io/rfcs/2000-const-generics.html) feature). +Interpolated trajectories are computed in parallel leading to a negligible impact in computational performance. To build a trajectory from a propagation event, call the [`for_duration_with_traj()`](https://rustdoc.nyxspace.com/nyx_space/propagators/struct.PropInstance.html#method.for_duration_with_traj) function on the propagator. Note that all of the event finding requires an interpolated trajectory. -??? info "Technical note" - Starting a propagator with the `for_duration_with_traj()` will return the final state of the propagation and the associated trajectory. This only works for states which also implement the [`InterpState` trait](https://docs.rs/nyx-space/latest/nyx_space/md/trajectory/trait.InterpState.html). +Trajectories in Nyx are a simple list of discrete states interpolated using the [Hermite interpolation](../../../anise/reference/mathspec/interpolation/hermite.md). They are akin to the NASA/SPICE Hermite Type 13 interpolator. - The segments in a trajectory are organized in a binary tree map allowing for blazing fast access to any segment of the trajectory. +Spacecraft trajectories can be transformed into another frame by transforming each of the individual states into the desired frame. -!!! warning "Limitations" - 1. To support generating the trajectories in parallell and storing them in a binary tree, it was necessary to choose an indexing method independent of the duration of the interpolation spline. This has been chosen to be the _centiseconds_ (10 milliseconds) since the start of the trajectory, the start being the very first state used to seed the trajectory. This means that if two splines last _less than 10 ms_, one will **overwrite** the other in the binary tree, meaning that part of the trajectory will **not** be queryable. - 1. Changing the frame of a trajectory follows the Shannon Nyquist sampling theorem. However, this will introduce _some imprecision_ nevertheless. Multiple conversions (e.g. from EME2000 to Moon J2000 back to EME2000) will accumulate errors quickly. - 1. Trajectories can currently only be generated with forward propagation. (Fixed in version 1.1.) +!!! warning + Transforming a trajectory into another frame _is different_ from propagating the original spacecraft in that other frame. A propagator will advance with a different step depending on the central body, refe to [the reference on propagator](../../../nyxspace/MathSpec/propagators.md). + + For example, propagating a spacecraft defined in the Earth Mean J2000 (EME2000) inertial frame when accounting for the point mass gravity of the Earth and the Moon and transforming that trajectory into a Moon J2000 centered frame will lead to a different trajectory that propagating that same inertial state defined in the Moon J2000 frame. -## Ephemeris +## Trajectory Ephemeris Integration Test -1. Propagate a LEO orbit with STM for 31 days in two body dynamics and request the propagator to generate the trajectory on the fly. -1. Compute the average SMA of that orbit by sampling the trajectory every 1 day (this ensures that the compiler does not remove the trajectory querying section from the test). -1. Check that the first and last state of the trajectory are strictly equal to the input propagator state and output propagator state. -1. Check that querying a trajectory one nanosecond after its end state returns an error. -1. Request a new propagation from the initial state but without generating the trajectory and read every intermediate state from that propagation. For each state, compute the radius and velocity differences between each of those states and the interpolated states at those exact times. -1. Ensure that the interpolation error in position is less than 1 micrometer and the error in velocity is less than 1 micrometer per second. -1. Convert the whole trajectory to its equivalent in the Moon J2000 frame (called `Luna` in Nyx). Convert that new trajectory back into the Earth Mean Equator J2000 frame (`EME2000`). -1. Sample that "double converted" trajectory every five minutes and compare each of those states with the initial interpolated trajectory. -1. Ensure that the interpolation error in position is less than 1 micrometer and the error in velocity is less than 1 micrometer per second. +Details of the `traj_ephem_forward` integration test. -All of the above is computed in 0.9 seconds on a standard desktop computer. +### Test Configuration -!!! note - If propagating with the STM, then the trajectory will also include the STM for that specific step forward. +- Initial epoch: 2021-01-01 12:00:00 UTC +- Initial state: LEO orbit in Earth J2000 frame +- Propagation duration: 31 days +- Dynamics: Two-body +- Integration method: variable step Runge-Kutta 8-9 -!!! check "Validation" - ```sh - $ RUST_LOG=info RUST_BACKTRACE=1 /usr/bin/time cargo test traj_ephem --release -- --nocapture - INFO nyx_space::propagators::propagator > Propagating for 30 days 23 h 59 min 60 s 0 ms -0.021928026217210607 ns until 2021-02-01T12:00:00 UTC - Average SMA: 7712.186 km Should be: 7712.186 - [tests/propagation/trajectory.rs:46] sum_sma / cnt - start_state.sma() = 0.000000245786395680625 - Ephem: Trajectory from 2021-01-01T12:00:00 UTC to 2021-02-01T12:00:00 UTC (30 days 23 h 59 min 60 s 0 ms -0.021928026217210607 ns, or 2678400.000 s) [4613 splines] - INFO nyx_space::propagators::propagator > Propagating for 30 days 23 h 59 min 60 s 0 ms -0.021928026217210607 ns until 2021-02-01T12:00:00 UTC - [traj_ephem] Maximum interpolation error: pos: 2.25e-8 m vel: 1.98e-11 m/s - INFO nyx_space::md::trajectory::traj > Converted trajectory from Earth J2000 to Moon J2000 in 189 ms - ephem_luna Trajectory from 2021-01-01T12:00:00 UTC to 2021-02-01T12:00:00 UTC (30 days 23 h 59 min 60 s 0 ms -0.021928026217210607 ns, or 2678400.000 s) [11203 splines] - INFO nyx_space::md::trajectory::traj > Converted trajectory from Moon J2000 to Earth J2000 in 448 ms - (...) +### Success Criteria - [traj_ephem] Maximum interpolation error after double conversion: pos: 1.27e-1 m vel: 4.41e-7 m/s - test propagation::trajectory::traj_ephem ... ok +| Metric | Requirement | +|--------|-------------| +| Position Error (Direct) | 0 m | +| Velocity Error (Direct) | 0 m/s | +| Position Error (Post frame conv.) | < 1 m | +| Velocity Error (Post frame conv.) | < 0.01 m/s | +| Temporal Precision | < 1 μs | - test result: ok. 1 passed; 0 failed; 1 ignored; 0 measured; 116 filtered out; finished in 0.94s - ``` +### Test Structure + +#### Trajectory Generation and Basic Validation + + - Builds a trajectory + - Validates conservation of semi-major axis + - Verifies initial and final states match propagation + - Ensures STM (State Transition Matrix) is properly unset + - Confirms interpolation bounds are respected + +#### Interpolation Accuracy + - Parallel truth generation + - Compares interpolated states against truth trajectory + - Validates position and velocity errors at exact timesteps + - Expected accuracy: machine precision for stored states -## Spacecraft trajectory +#### Persistence Testing -!!! check "Validation" + - Exports trajectory to Parquet format + - Includes eclipse event detection + - Reloads trajectory and verifies: + - State count consistency + - Temporal boundaries + - State vector equality + - Microsecond-level epoch precision + +#### Frame Transformation Validation + + - Double conversion test: Earth J2000 → Moon J2000 → Earth J2000 + - Samples every 5 minutes + - Validates: + - Temporal consistency + - Position error < 1 meter + - Velocity error < 0.01 m/s + + +!!! check "Run log" ```sh - $ RUST_LOG=info RUST_BACKTRACE=1 /usr/bin/time cargo test traj_spacecraft --release -- --nocapture + $ RUST_LOG=info RUST_BACKTRACE=1 /usr/bin/time cargo test traj_ephem --release -- --nocapture + INFO nyx_space::propagators::instance > Propagating for 31 days until 2021-02-01T12:00:00 UTC + INFO nyx_space::propagators::instance > ... done in 247 ms 822 μs 611 ns + [TIMING] 293 ms 488 μs 896 ns + Average SMA: 7712.186 km Should be: 7712.186 + Ephem: Trajectory in Earth J2000 (μ = 398600.435436096 km^3/s^2, eq. radius = 6378.14 km, polar radius = 6356.75 km, f = 0.0033536422844278) from 2021-01-01T12:00:00 UTC to 2021-02-01T12:00:00 UTC (31 days, or 2678400.000 s) [32289 states] + INFO nyx_space::propagators::instance > Propagating for 31 days until 2021-02-01T12:00:00 UTC + INFO nyx_space::propagators::instance > ... done in 274 ms 970 μs 367 ns + [traj_ephem] Maximum error on exact step: pos: 0.00e0 m vel: 0.00e0 m/s + INFO nyx_space::md::trajectory::traj > Exporting trajectory to parquet file... + INFO nyx_space::md::trajectory::traj > Serialized 32289 states from 2021-01-01T12:00:00 UTC to 2021-02-01T12:00:00 UTC + INFO nyx_space::md::trajectory::traj > Evaluating 1 event(s) + INFO nyx_space::md::trajectory::traj > Trajectory written to /home/chris/Workspace/nyx-space/nyx/output_data/ephem_forward-2024-10-27T06-36-19.parquet in 2 s 245 ms 95 μs 680 ns + INFO nyx_space::io::trajectory_data > File: /home/chris/Workspace/nyx-space/nyx/output_data/ephem_forward-2024-10-27T06-36-19.parquet + INFO nyx_space::io::trajectory_data > Created on: 2024-10-27T06:36:20.258980864 UTC + INFO nyx_space::io::trajectory_data > nyx-space License: AGPL 3.0 + INFO nyx_space::io::trajectory_data > Created by: Chris Rabotin (chris) on Linux + INFO nyx_space::io::trajectory_data > Purpose: Trajectory data + INFO nyx_space::io::trajectory_data > Generated by: Nyx v2.0.0-rc + INFO nyx_space::md::trajectory::sc_traj > Converted trajectory from Earth J2000 (μ = 398600.435436096 km^3/s^2, eq. radius = 6378.14 km, polar radius = 6356.75 km, f = 0.0033536422844278) to Moon J2000 in 300 ms: Trajectory in Moon J2000 (μ = 4902.800066163796 km^3/s^2, radius = 1737.4 km) from 2021-01-01T12:00:00 UTC to 2021-02-01T12:00:00 UTC (31 days, or 2678400.000 s) [32289 states] + ephem_luna Trajectory in Moon J2000 (μ = 4902.800066163796 km^3/s^2, radius = 1737.4 km) from 2021-01-01T12:00:00 UTC to 2021-02-01T12:00:00 UTC (31 days, or 2678400.000 s) [32289 states] + INFO nyx_space::md::trajectory::sc_traj > Converted trajectory from Moon J2000 (μ = 4902.800066163796 km^3/s^2, radius = 1737.4 km) to Earth J2000 (μ = 398600.435436096 km^3/s^2, eq. radius = 6378.14 km, polar radius = 6356.75 km, f = 0.0033536422844278) in 299 ms: Trajectory in Earth J2000 (μ = 398600.435436096 km^3/s^2, eq. radius = 6378.14 km, polar radius = 6356.75 km, f = 0.0033536422844278) from 2021-01-01T12:00:00 UTC to 2021-02-01T12:00:00 UTC (31 days, or 2678400.000 s) [32289 states] + Ephem back: Trajectory in Earth J2000 (μ = 398600.435436096 km^3/s^2, eq. radius = 6378.14 km, polar radius = 6356.75 km, f = 0.0033536422844278) from 2021-01-01T12:00:00 UTC to 2021-02-01T12:00:00 UTC (31 days, or 2678400.000 s) [32289 states] + 2.328666121016051e-11 + Eval: total mass = 0.000 kg @ [Earth J2000] 2021-01-01T12:05:00 UTC position = [-834.810409, -3848.218391, 6622.529651] km velocity = [5.519092, -4.261593, -1.778298] km/s Coast + Conv: total mass = 0.000 kg @ [Earth J2000] 2021-01-01T12:05:00 UTC position = [-834.810409, -3848.218391, 6622.529651] km velocity = [5.519092, -4.261593, -1.778298] km/s Coast 0.000 m + (...) - [traj_spacecraft] Maximum interpolation error: pos: 1.75e-7 m vel: 1.35e-6 m/s fuel: 1.14e-10 g full state: 3.06e-8 (no unit) - INFO nyx_space::md::trajectory::traj > Converted trajectory from Earth J2000 to Moon J2000 in 1 ms - INFO nyx_space::md::trajectory::traj > Converted trajectory from Moon J2000 to Earth J2000 in 3 ms - [traj_ephem] Maximum interpolation error after double conversion: pos: 5.76e-2 m vel: 1.52e-7 m/s - test traj_spacecraft ... ok - - test result: ok. 1 passed; 0 failed; 0 ignored; 0 measured; 2 filtered out; finished in 0.35s + [traj_ephem] Maximum interpolation error after double conversion: pos: 3.76e-8 m vel: 1.13e-9 m/s + test propagation::trajectory::traj_ephem_forward ... ok ``` -## Navigation trajectory - -Version 1.1 will allow generating navigation as well. \ No newline at end of file From 3390f2ae642e77fb1d5efc218a8565363682aa38 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sun, 27 Oct 2024 11:06:44 -0600 Subject: [PATCH 4/6] STM page updated --- .../MathSpec/celestial/coord_systems.md | 2 - docs/nyxspace/MathSpec/models/srp.md | 2 +- docs/nyxspace/MathSpec/optimization/stm.md | 63 +++++++++---------- mkdocs.yml | 2 +- 4 files changed, 33 insertions(+), 36 deletions(-) diff --git a/docs/nyxspace/MathSpec/celestial/coord_systems.md b/docs/nyxspace/MathSpec/celestial/coord_systems.md index 3a3c6da..dd0b456 100644 --- a/docs/nyxspace/MathSpec/celestial/coord_systems.md +++ b/docs/nyxspace/MathSpec/celestial/coord_systems.md @@ -24,8 +24,6 @@ $$ In practice, you can retrieve this rotation matrix by calling [`dcm_from_vnc_to_inertial`](https://docs.rs/anise/latest/anise/astro/orbit/type.Orbit.html#method.dcm_from_vnc_to_inertial) on an `Orbit` structure, or [`dcm3x3_from_vnc_to_inertial`](https://docs.rs/anise/latest/anise/astro/orbit/type.Orbit.html#method.dcm3x3_from_vnc_to_inertial) for the $3 \times 3$ DCM only. -Here are the two coordinate frame descriptions following the same style: - ### RIC Frame The **RIC frame** is a trajectory-centered coordinate system ideal for analyzing spacecraft motion relative to its position vector and orbital plane. It is particularly useful for relative navigation and rendezvous operations, where spacecraft position relative to a target is critical. diff --git a/docs/nyxspace/MathSpec/models/srp.md b/docs/nyxspace/MathSpec/models/srp.md index e03cc35..2207513 100644 --- a/docs/nyxspace/MathSpec/models/srp.md +++ b/docs/nyxspace/MathSpec/models/srp.md @@ -20,7 +20,7 @@ Note that $C_r$ and $\mathcal{A}$ are stored in the "context" passed to the EOM First, we compute the position of the Sun as seen from the spacecraft, and its unit vector, respectively $\mathbf{r_\odot}$ and $\mathbf{\hat r_\odot}$. Then, we compute the shadowing factor, $k$, using the [eclipse model](../celestial/eclipse.md). -Compute the norm of Sun vector in AU, $||\mathbf{r_\odot}||_{\text{AU}}$ by dividing the $\mathbf{r_\odot}$ vector by 1 AU. +Compute the norm of Sun vector in AU, $ ||\mathbf{r_\odot}||_{\text{AU}}$ by dividing the $\mathbf{r_\odot} $ vector by 1 AU. Compute the flux pressure as follows: diff --git a/docs/nyxspace/MathSpec/optimization/stm.md b/docs/nyxspace/MathSpec/optimization/stm.md index b986451..4e4751f 100644 --- a/docs/nyxspace/MathSpec/optimization/stm.md +++ b/docs/nyxspace/MathSpec/optimization/stm.md @@ -1,48 +1,47 @@ -# State Transition Matrix computation +The State Transition Matrix (STM, $\boldsymbol \Phi$) is a linearization tool for dynamical systems: it is a local approximation of the dynamics of the system. -Nyx uses dual number theory (i.e. automatic differentiation) to compute all state transition matrices (or Jacobian matrices) to machine precision without relying on finite differencing or error-prone manual derivations. All of the computations are handled by the [hyperdual](https://gitlab.com/chrisrabotin/hyperdual) Rust package, which I co-authored. +In the field of astrodynamics, it is an approximation of the forces applied to the spacecraft over a short period of time. The STM is a necessary tool in orbit determination and in the circular restricted three-body problem (CRTBP), and may also be used in trajectory optimization. -The following is an extract of the AAS-19-716 conference paper called _Application of Dual Number Theory to Statistical Orbital Determination_, by C. Rabotin [^1]. +For example, in orbit determination, the STM is computed from a propagated spacecraft state. It is considered a reasonable approximation of the forces acting on a spacecraft until the next time step (typically the next tracking measurement) through the time-propagation of the whole matrix. -First, the mathematical notions of dual numbers and hyper-dual spaces is recalled, along with their error-free auto-differentiation properties. Subsequently, hyper-dual spaces are elegantly applied to statistical orbital determination problems, both for the computation of linearized dynamics, and for the sensitivity matrix of a Kalman filter. Further, a comparison of the execution speeds of the same statistical orbital determination problem using the hyper-dual space formulation and the traditional manual derivation formulation is provided. +## System Dynamics -Please refer to [Dual Numbers](../appendix/dual_numbers.md) for a primer on dual number theory. +Consider an orbital state $\boldsymbol{X}$ in a Cartesian frame, comprising position $\boldsymbol{r} = [x,y,z]$ and velocity $\boldsymbol{v} = [\dot{x}, \dot{y}, \dot{z}]$: -## Applying hyper-dual space to statistical orbital determination using Kalman filtering -### The state transition matrix -The state transition matrix (STM, $\boldsymbol \Phi$) is a linearization procedure of a dynamical system. In the field of astrodynamics, it is an approximation of the dynamics over a short period of time. In the case of statistical orbital determination, the STM is computed around a reference point of an orbit and propagated forward in time until the next propagation time step or until the next spacecraft observation. The STM is computed via the partials matrix $\boldsymbol A$ of the state $\boldsymbol X$ at a time $t$, as shown earlier. +$$\boldsymbol{X} = \begin{bmatrix} x \\ y \\ z \\ \dot{x} \\ \dot{y} \\ \dot{z} \end{bmatrix}$$ -\begin{equation} -\label{grad_stm} -\boldsymbol A(t) = \frac{d \boldsymbol X(t)}{d\boldsymbol X_0} -\end{equation} +The time derivative $\boldsymbol{\dot{X}}$ represents system dynamics: -If one only estimates the spacecraft state, then the partials matrix corresponds to the gradient of the orbital dynamics applied to the spacecraft. The state $\boldsymbol X$ of a spacecraft may be defined as the vector of its position and velocity in a Cartesian frame. The equation below shows the relation between the partials matrix and the STM, where $t$ and $t_0$ correspond to two different times such that $t>t_0$. +$$\boldsymbol{\dot{X}} = F(\boldsymbol{X}, t)$$ -\begin{equation} -\label{stm_def} -\frac{d {\boldsymbol \Phi}(t, t_0)}{dt} = \boldsymbol A(t)\boldsymbol \Phi(t, t_0) -\end{equation} +For two-body problems, acceleration $\boldsymbol{\dot{v}}$ is $-\frac{\mu}{r^3}\boldsymbol{X}$. More complex dynamics extend $F$ to include additional state components. For example, when enabling [Solar Radiation Pressure (SRP)](../models/srp.md) estimation in Nyx, the state vector is expanded to include the coefficient of reflectivity $C_r$. -The $\boldsymbol A$ matrix contains the partial derivatives of the accelerations ($a_x,~a_y,~a_z$) with respect to each of the components of the position of the spacecraft, cf. the next equation. +### Computation -\begin{equation} -\label{tb_partials} -\boldsymbol A = \begin{bmatrix} - 0 & 0 & 0 & 1 & 0 & 0 \\ - 0 & 0 & 0 & 0 & 1 & 0 \\ - 0 & 0 & 0 & 0 & 0 & 1 \\ - \frac{\partial a_x}{\partial x} & \frac{\partial a_y}{\partial x} & \frac{\partial a_z}{\partial x} & 0 & 0 & 0 \\ - \frac{\partial a_x}{\partial y} & \frac{\partial a_y}{\partial y} & \frac{\partial a_z}{\partial y} & 0 & 0 & 0 \\ - \frac{\partial a_x}{\partial z} & \frac{\partial a_y}{\partial z} & \frac{\partial a_z}{\partial z} & 0 & 0 & 0 \\ -\end{bmatrix} -\end{equation} +The STM is computed via the partials matrix (or _Jacobian_) $\boldsymbol A$ of the state $\boldsymbol X$ at a time $t$. A first order Taylor series expansion is used to compute this Jacobian. -In the case of two-body dynamics, the $\boldsymbol A$ matrix is relatively simple to compute. Deriving this partials matrix for higher fidelity dynamics or for a larger state is more complicated. In practice, this may lead orbital estimation software to ignore part of these dynamics. Navigators may then account for smaller perturbations by inflating the covariance matrix with stochastic noise. +$$ \boldsymbol{\dot{X}} = \boldsymbol{\dot{X}^*} + \frac{\partial F (\boldsymbol{X^*})}{\partial \boldsymbol{\dot{X}}} \cdot \left( \boldsymbol{\dot{X}} - \boldsymbol{\dot{X}^*} \right) $$ -Dual numbers, however, provide the partials matrix as part of the computation of the equations of motion (EOM) themselves, as long as these EOMs describe the movement of all of the state variables to be estimated. In practice, this requires defining a hyper-dual space whose size is equal to the number of variables in the state to be estimated. For example, if estimating the Cartesian state a spacecraft and a maneuver magnitude at a reference point during an orbital determination arc, building a seven-dimensional dual space will return the result of the EOMs and the partials matrix computed at that point. An open source example of the building of such hyper-dual space has been implemented as a test case in \textit{dual\_num}, a thorough dual number library in a programming language called Rust. +where $\boldsymbol{X^*}$ is a reference state vector, e.g. the state at $t_0$. +In practice, this Jacobian is the matrix of partial derivatives of the acceleration of the system, forming the $A$ matrix, where $\boldsymbol{a_{x,y,z}}$ is the acceleration in the X, Y, and Z axes: -[^1]: That's me. `\o/` +$$\boldsymbol A = \frac{d \boldsymbol X(t)}{d\boldsymbol X_0} = \begin{bmatrix} 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ \frac{\partial a_x}{\partial x} & \frac{\partial a_y}{\partial x} & \frac{\partial a_z}{\partial x} & 0 & 0 & 0 \\ \frac{\partial a_x}{\partial y} & \frac{\partial a_y}{\partial y} & \frac{\partial a_z}{\partial y} & 0 & 0 & 0 \\ \frac{\partial a_x}{\partial z} & \frac{\partial a_y}{\partial z} & \frac{\partial a_z}{\partial z} & 0 & 0 & 0 \end{bmatrix}$$ + +## Dual Numbers and Hyper-Dual Spaces + +Nyx computes STMs using dual number theory because it avoids manual derivation errors via its error-free auto-differentiation properties. This is implemented via the [hyperdual](https://github.com/christopherrabotin/hyperdual) Rust crate. Refer to the page on [Dual Number theory](../appendix/dual_numbers.md) for a primer on the topic. + +### Application in Kalman Filtering and Orbit determination + +In statistical orbital determination, the STM is computed around an orbit reference point and propagated forward. The partials matrix $\boldsymbol A(t)$ relates to the STM as: + +$$\frac{d {\boldsymbol \Phi}(t, t_0)}{dt} = \boldsymbol A(t)\boldsymbol \Phi(t, t_0)$$ + +Refer to the [Kalman filtering](../orbit_determination/kalman.md) page for details on how the state transition matrix is used in orbit determination in Nyx. + +Dual numbers provide these derivatives directly from equations of motion, enabling efficient computation even in complex systems. For details, refer to the AAS-19-716 conference paper _Application of Dual Number Theory to Statistical Orbital Determination_, Rabotin. + +[^1]: For detailed STM discussion, refer to section 1.2.5 and 4.2.1. "Statistical Orbit Determination" by Tapley et al., Elsevier 2004. --8<-- "includes/Abbreviations.md" \ No newline at end of file diff --git a/mkdocs.yml b/mkdocs.yml index b14f375..5fd2d13 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -73,7 +73,7 @@ nav: - Overview: nyxspace/MathSpec/monte_carlo/overview.md - Error modeling: nyxspace/MathSpec/monte_carlo/error_models.md - Optimization: - - State transition matrix: nyxspace/MathSpec/optimization/stm.md + - State transition matrix (STM): nyxspace/MathSpec/optimization/stm.md - Targeting: nyxspace/MathSpec/optimization/targeting.md - Propagators: nyxspace/MathSpec/propagators.md - Appendices: From de12cacdaa2a6c613414f211daf1b48beddc0d56 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sun, 27 Oct 2024 11:17:14 -0600 Subject: [PATCH 5/6] Fix SRP page --- docs/nyxspace/MathSpec/models/srp.md | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/docs/nyxspace/MathSpec/models/srp.md b/docs/nyxspace/MathSpec/models/srp.md index 2207513..8a43552 100644 --- a/docs/nyxspace/MathSpec/models/srp.md +++ b/docs/nyxspace/MathSpec/models/srp.md @@ -1,5 +1,3 @@ -# Solar radiation pressure - SRP is implemented as a `ForceModel`, meaning that the implementation computes a force from which the instantaneous mass of the spacecraft is subsequently divided to apply an acceleration. Therefore, as opposed to most formulations in the literature, we will use a force formulation. Currently, Nyx only supports spherical SRP model. @@ -20,18 +18,19 @@ Note that $C_r$ and $\mathcal{A}$ are stored in the "context" passed to the EOM First, we compute the position of the Sun as seen from the spacecraft, and its unit vector, respectively $\mathbf{r_\odot}$ and $\mathbf{\hat r_\odot}$. Then, we compute the shadowing factor, $k$, using the [eclipse model](../celestial/eclipse.md). -Compute the norm of Sun vector in AU, $ ||\mathbf{r_\odot}||_{\text{AU}}$ by dividing the $\mathbf{r_\odot} $ vector by 1 AU. +Compute the norm of Sun vector in AU, $||\mathbf{r_\odot}||$ by dividing the $\mathbf{r_\odot}$ vector by 1 AU. -Compute the flux pressure as follows: +Compute the flux pressure ($\Phi$) as follows: -$$\Phi_{\text{SRP}} = \frac{k\phi}{c} \left(\frac{1.0}{||\mathbf{r_\odot}||_{\text{AU}}} \right)^2$$ +$$\Phi = \frac{k\phi}{c} \left(\frac{1.0}{||\mathbf{r_\odot}||} \right)^2$$ Finally, return the SRP force [^2]: -$$ \mathbf{F}_{\text{SRP}} = C_r \mathcal{A} \Phi_{\text{SRP}} \mathbf{\hat r_\odot}$$ +$$ \mathbf{F} = C_r \mathcal{A} \Phi \mathbf{\hat r_\odot}$$ + +## Estimation -!!! note - Although the above derivation mentions the Sun, Nyx trivially supports any other light source regardless of the integration frame. +Nyx can estimate the coefficient of reflectivity $C_r$ in orbit determination scenarios, as shown in [the Lunar Reconnaissance Orbiter](../../showcase/04_lro_od/index.md) example. ??? check "Validation" Nyx has three validation scenarios for the SRP computation to ensure that we test full illumination (`srp_earth_full_vis`), long penumbra passages (`srp_earth_meo_ecc_inc`), and very short penumbra passages (`srp_earth_penumbra`). In all of the test cases, we propagate a spacecraft for 24 **days** to ensure that a high amount of error can accumulate if the modeling is incorrect. The worst absolute position error is _high_ compared to GMAT: 287 meters. From d7ef619e92cec5566c1096679bf3bf37078cec9c Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sun, 27 Oct 2024 11:26:59 -0600 Subject: [PATCH 6/6] Good to go --- .../reference/mathspec/interpolation/chebyshev.md | 10 ++++++++++ docs/anise/reference/mathspec/interpolation/index.md | 6 +++++- docs/nyxspace/MathSpec/celestial/interpolation.md | 6 +++--- 3 files changed, 18 insertions(+), 4 deletions(-) diff --git a/docs/anise/reference/mathspec/interpolation/chebyshev.md b/docs/anise/reference/mathspec/interpolation/chebyshev.md index 94bec35..bbebe8b 100644 --- a/docs/anise/reference/mathspec/interpolation/chebyshev.md +++ b/docs/anise/reference/mathspec/interpolation/chebyshev.md @@ -52,6 +52,16 @@ For $j = N \to 2$: $$ f'(x) = \frac{w_0 + xdw_0 - dw_1}{r_s} $$ +## Output + +The algorithm returns a tuple $(f, f')$ where: + +- $f$ is the interpolated function value at $x$ +- $f'$ is the interpolated derivative at $x$ + +!!! note + A variation of this algorithm is used for Type 3 SPK/BSP SPICE files where the time derivative is not computed. + ## Error Handling The implementation checks for: diff --git a/docs/anise/reference/mathspec/interpolation/index.md b/docs/anise/reference/mathspec/interpolation/index.md index b077e6f..ada894d 100644 --- a/docs/anise/reference/mathspec/interpolation/index.md +++ b/docs/anise/reference/mathspec/interpolation/index.md @@ -1 +1,5 @@ -ANISE, like SPICE, heavily relies on interpolation for ephemeris and attitude trajectories. This section details each of the algorithms used in ANISE. \ No newline at end of file +ANISE, like SPICE, heavily relies on interpolation for ephemeris and attitude trajectories. This section details each of the algorithms used in ANISE. + +- [Chebyshev](./chebyshev.md) +- [Hermite](./hermite.md) +- [Lagrange](./lagrange.md) \ No newline at end of file diff --git a/docs/nyxspace/MathSpec/celestial/interpolation.md b/docs/nyxspace/MathSpec/celestial/interpolation.md index b0a4571..7812f8c 100644 --- a/docs/nyxspace/MathSpec/celestial/interpolation.md +++ b/docs/nyxspace/MathSpec/celestial/interpolation.md @@ -4,14 +4,14 @@ For details on interpolation methods for trajectories, including [Chebyshev](../ Interpolated trajectories are computed in parallel leading to a negligible impact in computational performance. To build a trajectory from a propagation event, call the [`for_duration_with_traj()`](https://rustdoc.nyxspace.com/nyx_space/propagators/struct.PropInstance.html#method.for_duration_with_traj) function on the propagator. Note that all of the event finding requires an interpolated trajectory. -Trajectories in Nyx are a simple list of discrete states interpolated using the [Hermite interpolation](../../../anise/reference/mathspec/interpolation/hermite.md). They are akin to the NASA/SPICE Hermite Type 13 interpolator. +Trajectories in Nyx are a simple list of discrete states interpolated using the [Hermite interpolator](../../../anise/reference/mathspec/interpolation/hermite.md). They are akin to the NASA/SPICE Hermite Type 13 interpolator. Spacecraft trajectories can be transformed into another frame by transforming each of the individual states into the desired frame. !!! warning - Transforming a trajectory into another frame _is different_ from propagating the original spacecraft in that other frame. A propagator will advance with a different step depending on the central body, refe to [the reference on propagator](../../../nyxspace/MathSpec/propagators.md). + Transforming a trajectory into another frame _is different_ from propagating the original spacecraft in that other frame. A propagator will advance with a different step depending on the central body, refer to [the reference on propagator](../../../nyxspace/MathSpec/propagators.md). - For example, propagating a spacecraft defined in the Earth Mean J2000 (EME2000) inertial frame when accounting for the point mass gravity of the Earth and the Moon and transforming that trajectory into a Moon J2000 centered frame will lead to a different trajectory that propagating that same inertial state defined in the Moon J2000 frame. + For example, propagating a spacecraft defined in the Earth J2000 inertial frame when accounting for the gravitational forces of the Earth and the Moon, and transforming that trajectory into a Moon J2000 centered frame will lead to a different trajectory that propagating that same inertial state defined in the Moon J2000 frame. ## Trajectory Ephemeris Integration Test