Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 6 additions & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
85 changes: 85 additions & 0 deletions docs/anise/reference/mathspec/interpolation/chebyshev.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,85 @@
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} $$

## 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:

- 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.
77 changes: 77 additions & 0 deletions docs/anise/reference/mathspec/interpolation/hermite.md
Original file line number Diff line number Diff line change
@@ -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.
5 changes: 5 additions & 0 deletions docs/anise/reference/mathspec/interpolation/index.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
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)
71 changes: 71 additions & 0 deletions docs/anise/reference/mathspec/interpolation/lagrange.md
Original file line number Diff line number Diff line change
@@ -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.
4 changes: 2 additions & 2 deletions docs/blog/posts/2024.10.18-hifitime-4.0.0.md
Original file line number Diff line number Diff line change
Expand Up @@ -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 }

Expand Down
2 changes: 1 addition & 1 deletion docs/hifitime/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down
Loading
Loading