James Webb Space Telescope covariance Monte Carlo¶
Spaceflight dynamics requires a lot of statistical analysis to convince ourselves that the results we'll present to all other teams are correct.
In this example, you'll learn how to sample the orbital state from SPICE BSP file (the industry standard for ephemeris information) and build a spacecraft structure around that orbit. Then, we'll build an uncertainty in the position and velocity of that spacecraft, and propagate it into the future.
To run this example, just execute:
Building in release mode will make the computation significantly faster. Specifying RUST_LOG=info will allow you to see all of the information messages happening in ANISE and Nyx throughout the execution of the program.
Objective¶
We aim to compare two different statistical methods in order to demonstrate two advanced features of Nyx.
First, we'll use a covariance mapping approach, whereby the covariance at some time t0 is propagated to a future time t1 along with the spacecraft trajectory. This is commonly used for the time update or prediction step of an orbit determination process.
Then, we'll use the Monte Carlo framework of Nyx to propagate the initial spacecraft state after dispersions using all threads of the computer. Multi-threaded propagation is not common in other astrodynamics software.
Finally, we'll check that the 3-sigma (i.e. 99.7%) covariance bounds of the covariance mapping approach matches the Monte Carlo results in terms of uncertainty in the state vector and in the Keplerian orbital elements.
Example run¶
Nyx is blazing fast. The covariance mapping and the Monte Carlo runs of 5000 runs is executed in less than one minute. Then it takes about 41 seconds to export all of the trajectory data into a ~246 MB parquet file.
Compiling nyx-space v2.0.0-rc (/home/crabotin/Workspace/nyx-space/nyx)
Finished `release` profile [optimized] target(s) in 3.52s
Running `target/release/examples/02_jwst`
INFO anise::almanac::metaload::metafile > Saved https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/jwst_rec.bsp to /home/crabotin/.local/share/nyx-space/anise/jwst_rec.bsp (CRC32 = a8460057)
INFO anise::almanac::metaload::metafile > Using cached /home/crabotin/.local/share/nyx-space/anise/de440s.bsp
INFO anise::almanac::metaload::metafile > Using cached /home/crabotin/.local/share/nyx-space/anise/pck11.pca
INFO anise::almanac::metaload::metafile > Discarding cached /home/crabotin/.local/share/nyx-space/anise/moon_fk.epa - CRC32 differ (got 194230817, config expected 292928914)
INFO anise::almanac::metaload::metafile > Saved http://public-data.nyxspace.com/anise/v0.4/moon_fk.epa to /home/crabotin/.local/share/nyx-space/anise/moon_fk.epa (CRC32 = b93ba21)
INFO anise::almanac::metaload::metafile > Discarding cached /home/crabotin/.local/share/nyx-space/anise/moon_pa_de440_200625.bpc - CRC32 differ (got 3454388861, config expected 1817759242)
INFO anise::almanac::metaload::metafile > Saved http://public-data.nyxspace.com/anise/moon_pa_de440_200625.bpc to /home/crabotin/.local/share/nyx-space/anise/moon_pa_de440_200625.bpc (CRC32 = cde5ca7d)
INFO anise::almanac::metaload::metafile > Saved https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/earth_latest_high_prec.bpc to /home/crabotin/.local/share/nyx-space/anise/earth_latest_high_prec.bpc (CRC32 = a74b6afd)
INFO anise::almanac > Loading almanac from /home/crabotin/.local/share/nyx-space/anise/de440s.bsp
INFO anise::almanac > Loading as DAF/SPK
INFO anise::almanac > Loading almanac from /home/crabotin/.local/share/nyx-space/anise/pck11.pca
INFO anise::almanac > Loading almanac from /home/crabotin/.local/share/nyx-space/anise/moon_fk.epa
INFO anise::almanac > Loading almanac from /home/crabotin/.local/share/nyx-space/anise/moon_pa_de440_200625.bpc
INFO anise::almanac > Loading as DAF/PCK
INFO anise::almanac > Loading almanac from /home/crabotin/.local/share/nyx-space/anise/earth_latest_high_prec.bpc
INFO anise::almanac > Loading as DAF/PCK
INFO anise::almanac > Loading almanac from /home/crabotin/.local/share/nyx-space/anise/jwst_rec.bsp
INFO anise::almanac > Loading as DAF/SPK
JWST defined from 2024-04-13T22:40:09.185630849 ET to 2024-07-08T00:01:09.183912447 ET
[Earth J2000] 2024-07-08T00:01:09.183912447 ET sma = 876856.027357 km ecc = 0.985029 inc = 52.340344 deg raan = 310.126394 deg aop = 130.798599 deg ta = 180.103906 deg
total mass = 6200.000 kg @ [Earth J2000] 2024-07-08T00:01:09.183912447 ET position = [119901.070276, -1389299.665421, -1041369.150539] km velocity = [0.045956, -0.013168, 0.034535] km/s Coast
RIC σ_x = 0.5 km σ_y = 0.3 km σ_z = 1.5 km
RIC σ_vx = 0.0001 km/s σ_vy = 0.0006 km/s σ_vz = 0.003 km/s
σ_cr = 0 σ_cd = 0 σ_mass = 0 kg
INFO nyx_space::od::process > Mapping covariance for 6 days 12 h with 1 min step
INFO nyx_space::od::process::export > Exporting orbit determination result to parquet file...
INFO nyx_space::od::process::export > Serialized 9360 estimates and residuals
INFO nyx_space::od::process::export > Orbit determination results written to ./02_jwst_covar_map.parquet in 734 ms 125 μs 312 ns
INFO nyx_space::mc::montecarlo > Propagated 5000 states in 58 s 792 ms 823 μs 694 ns
INFO nyx_space::mc::results > Exporting Monte Carlo results to parquet file...
INFO nyx_space::mc::results > Serialized 1065006 states from 2024-06-24T00:01:09.184303103 ET to 2024-06-30T12:01:09.184303103 ET
INFO nyx_space::mc::results > Evaluating 2 event(s)
INFO nyx_space::mc::results > Trajectory written to 02_jwst_monte_carlo.parquet in 41 s 603 ms 696 μs 128 ns
Analysis¶
Overall, we can confirm that the 3-sigma covariance is a good approximation of the uncertainty. Notably, all of the dispersed trajectories fit (nearly) within the 3-σ bounds in the state space and in the Keplerian orbital element space. The rotation from the Keplerian orbital element space into the Cartesian space is done by building the Jacobian from the desired Keplerian elements into the Cartesian space using hyperdual numbers (also called "automatic differentiation" in the machine learning world). This ensure that we don't lose any precision in the state space change, and that's shown to be the case in these plots.
State uncertainties¶
As expected from any orbit determination software, Nyx can output uncertainties in the state vector in the integration frame and in the RIC frame. Note: these plots look pretty linear, but that's because we're running a pure prediction filter and JWST is in a stable halo orbit.
Keplerian uncertainties¶
A few tools try to provide Keplerian uncertainties, but often fail to do so correctly (I'm looking at you, ODTK). Nyx rotates the covariance from its Cartesian form into the Keplerian state space by computing the partial derivatives of each requested parameter with respect to the nominal state. This computation is flawless because it uses automatic differentiation (via hyperdual numbers). As such, the OD export also includes all of the state computations supported in Nyx, including uncommon ones like the uncertainties in the energy of the orbit or in the true anomaly.
The following columns are also provided as 1-sigma in the OD dataframe: - Sigma aol (deg) - Sigma c3 (km^2/s^2) - Sigma declin (deg) - Sigma ea (deg) - Sigma fpa (deg) - Sigma hmag (km) - Sigma hx (km) - Sigma hy (km) - Sigma hz (km) - Sigma ma (deg) - Sigma right_asc (deg) - Sigma rmag (km) - Sigma semi_parameter (km) - Sigma semi_minor (km) - Sigma tlong (deg) - Sigma vmag (km/s) - Sigma Cr (Earth J2000) (unitless) - Sigma Cd (Earth J2000) (unitless) - Sigma Mass (Earth J2000) (kg)
