.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_stm.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_plot_stm.py: State Transition Matrix ======================= The State Transition Matrix (STM) describes how a small change in an object's initial state (position and velocity) maps to a change in its state at a later time: .. math:: \delta \mathbf{x}(t_f) \approx \Phi(t_f, t_0) \, \delta \mathbf{x}(t_0) This is the foundation of linear orbit determination: rather than re-integrating the orbit for every trial initial condition, a single STM integration gives a first-order approximation valid for small perturbations. kete computes the STM using the Radau 15th-order integrator with full N-body physics (all planets, GR, J2 oblateness). For objects with non-gravitational forces (comets, dust), extra columns give the partial derivatives of the final state with respect to the non-grav parameters, enabling simultaneous fitting of the orbit and the force model. This example shows: 1. Computing the 6x6 STM for a simple asteroid orbit. 2. Using the STM to propagate an orbital covariance matrix. 3. Visualizing the 1-sigma position uncertainty ellipse over time. .. GENERATED FROM PYTHON SOURCE LINES 28-34 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import kete .. GENERATED FROM PYTHON SOURCE LINES 35-41 Define an initial state and covariance --------------------------------------- We use Ceres as a convenient real object with a well-known orbit, then assign a synthetic diagonal covariance to give it ~10 km position and ~1 m/s velocity uncertainty (expressed in AU and AU/day). .. GENERATED FROM PYTHON SOURCE LINES 41-51 .. code-block:: Python jd_start = kete.Time.j2000().jd state = kete.spice.get_state("Ceres", jd_start) # 1-sigma uncertainties sigma_pos_au = 10 / 1.496e8 # 10 km in AU sigma_vel_auday = 1 / 1.731e6 # 1 m/s in AU/day cov0 = np.diag([sigma_pos_au**2] * 3 + [sigma_vel_auday**2] * 3) .. GENERATED FROM PYTHON SOURCE LINES 52-62 Propagate covariance over 1 year --------------------------------- Rather than re-propagating from the initial state at each time step, we chain each integration: the state and covariance from step N become the inputs to step N+1. The covariance update at each step is: .. math:: P(t_{k+1}) = \Phi_k \, P(t_k) \, \Phi_k^T .. GENERATED FROM PYTHON SOURCE LINES 62-79 .. code-block:: Python n_steps = 60 step_days = 365.0 / n_steps cur_state = state cur_cov = cov0 pos_sigma_au = [np.sqrt(np.trace(cur_cov[:3, :3]))] time_steps = [0.0] for k in range(n_steps): jd_next = cur_state.jd + step_days cur_state, stm = kete.state_transition.compute_stm(cur_state, jd_next) cur_cov = stm @ cur_cov @ stm.T pos_sigma_au.append(np.sqrt(np.trace(cur_cov[:3, :3]))) time_steps.append((k + 1) * step_days) pos_sigma_km = np.array(pos_sigma_au) * 1.496e8 .. GENERATED FROM PYTHON SOURCE LINES 80-85 Inspect the STM at 30 days --------------------------- ``compute_stm`` returns the propagated state and the full 6x6 sensitivity matrix. The STM must be symplectic (det ~= 1) for conservative dynamics. .. GENERATED FROM PYTHON SOURCE LINES 85-92 .. code-block:: Python jd_30 = jd_start + 30 final_state, stm = kete.state_transition.compute_stm(state, jd_30) print(f"STM determinant at 30 days: {np.linalg.det(stm):.8f} (should be ~1.0)") print(f"Final state: {final_state}") .. rst-class:: sphx-glr-script-out .. code-block:: none STM determinant at 30 days: 1.00000000 (should be ~1.0) Final state: State(desig="ceres", jd=2451575.0, pos=[-2.467541144091866, 0.47238990295848526, 0.4692916232900171], vel=[-0.0022900074001403485, -0.01092475456908372, 8.488499757221876e-5], frame=Ecliptic, center="sun") .. GENERATED FROM PYTHON SOURCE LINES 93-100 Sensitivity matrix with a non-gravitational model -------------------------------------------------- For a comet-like object, we can fit A1 (radial), A2 (transverse), and A3 (normal) non-grav parameters simultaneously with the orbit. The STM returns a 6x9 matrix: the first 6 columns are the standard STM, the last 3 are d(final state)/d(A1, A2, A3). .. GENERATED FROM PYTHON SOURCE LINES 100-112 .. code-block:: Python ng_model = kete.propagation.NonGravModel.new_comet( a1=1e-9, # AU/day^2 (typical weak cometary force) a2=3e-10, a3=0.0, ) _, sens = kete.state_transition.compute_stm(state, jd_30, non_grav=ng_model) print(f"\nSensitivity matrix shape with JplComet model: {sens.shape} (expected: 6x9)") print("Position sensitivity to A1 at 30 days (AU per AU/day^2):", sens[:3, 6]) .. rst-class:: sphx-glr-script-out .. code-block:: none Sensitivity matrix shape with JplComet model: (6, 9) (expected: 6x9) Position sensitivity to A1 at 30 days (AU per AU/day^2): [-6.37638895 1.18048605 1.8537434 ] .. GENERATED FROM PYTHON SOURCE LINES 113-115 Plot position uncertainty over time ------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 115-125 .. code-block:: Python fig, ax = plt.subplots(figsize=(8, 4)) ax.plot(time_steps, pos_sigma_km) ax.set_xlabel("Days from J2000") ax.set_ylabel("RMS position 1-sigma (km)") ax.set_title("Ceres orbit uncertainty propagation (full N-body STM)") ax.set_yscale("log") ax.grid(True, which="both", alpha=0.3) plt.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_stm_001.png :alt: Ceres orbit uncertainty propagation (full N-body STM) :srcset: /auto_examples/images/sphx_glr_plot_stm_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.258 seconds) .. _sphx_glr_download_auto_examples_plot_stm.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_stm.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_stm.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_stm.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_