Skip to content

Variational Equations and State Transition Matrix

Variational equations enable propagating not just the statye dynamics, but also how small perturbations would affect the state over time. They help relate how changes in state at one time map to changes at a later time which is critical for orbit determination and control.

What are Variational Equations?

For a dynamical system \(\dot{\mathbf{x}} = \mathbf{f}(t, \mathbf{x})\), variational equations describe how small deviations from the nominal trajectory evolve. Consider two nearby initial conditions:

  • Nominal: \(\mathbf{x}_0\)
  • Perturbed: \(\mathbf{x}_0 + \delta\mathbf{x}_0\)

The difference in trajectories can be approximated by the State Transition Matrix (STM) \(\Phi(t, t_0)\):

\[\delta\mathbf{x}(t) \approx \Phi(t, t_0) \cdot \delta\mathbf{x}_0\]

This relationship is exact for linear systems and accurate for nonlinear systems when \(||\delta\mathbf{x}_0||\) is small.

The State Transition Matrix

The State Transition Matrix (STM) satisfies the matrix differential equation:

\[\dot{\Phi}(t, t_0) = \mathbf{J}(t, \mathbf{x}(t)) \cdot \Phi(t, t_0)\]

where \(\mathbf{J}\) is the Jacobian matrix of the dynamics:

\[\mathbf{J}_{ij} = \frac{\partial f_i}{\partial x_j}\]

The STM has a few key properites. First, the intial condition of the STM is always the indentity matrix. That is:

\(\Phi(t_0, t_0) = \mathbf{I}\) (identity matrix)

For linear systems: \(\Phi(t, t_0)\) is the matrix exponential \(e^{\mathbf{A}(t-t_0)}\).

STM Propagation in Brahe

Brahe integrators can propagate the state and STM simultaneously using step_with_varmat(). What happens under the hood is:

  1. The integrator advances both the state (\(\mathbf{x}\)) and the STM (\(\Phi\)) using the same time step
  2. At each stage of the Runge-Kutta method, the Jacobian is evaluated at the current state
  3. The variational equations \(\dot{\Phi} = \mathbf{J} \cdot \Phi\) are integrated alongside the state equations
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
import brahe as bh
import numpy as np

# Initialize EOP
bh.initialize_eop()


# Define two-body dynamics
def dynamics(t, state):
    mu = bh.GM_EARTH
    r = state[0:3]
    v = state[3:6]
    r_norm = np.linalg.norm(r)
    a = -mu / r_norm**3 * r
    return np.concatenate([v, a])


# Create numerical Jacobian for variational equations
jacobian = bh.NumericalJacobian.central(dynamics).with_adaptive(
    scale_factor=1.0, min_value=1e-6
)

# Initial orbit (LEO)
r0 = np.array([bh.R_EARTH + 600e3, 0.0, 0.0])
v0 = np.array([0.0, 7.5e3, 0.0])
state0 = np.concatenate([r0, v0])

# Initial state transition matrix (identity)
phi0 = np.eye(6)

print("Integration with State Transition Matrix")
print(f"Initial orbit: {r0[0] / 1e3:.1f} km altitude")

# Create integrator with Jacobian
config = bh.IntegratorConfig.adaptive(abs_tol=1e-12, rel_tol=1e-11)
integrator = bh.DP54Integrator(6, dynamics, jacobian=jacobian, config=config)

# Propagate for one orbit period
t = 0.0
state = state0.copy()
phi = phi0.copy()
dt = 60.0

# Approximate orbital period
period = bh.orbital_period(np.linalg.norm(r0))

print("Time   Position STM[0,0]  Velocity STM[3,3]  Det(STM)")
print("-" * 60)

steps = 0
while t < period:
    # Propagate state and STM together (adaptive integrator returns 5-tuple)
    new_state, new_phi, dt_used, error_est, dt_next = integrator.step_with_varmat(
        t, state, phi, min(dt, period - t)
    )

    t += dt_used
    state = new_state
    phi = new_phi
    dt = dt_next
    steps += 1

    # Print progress
    if steps % 20 == 1:
        det_phi = np.linalg.det(phi)
        print(
            f"{t:6.0f}s    {phi[0, 0]:8.4f}      {phi[3, 3]:8.4f}        {det_phi:8.4f}"
        )

print(f"\nPropagation complete! ({steps} steps)")

# Example: Map initial position uncertainty to final uncertainty
print("\nExample: Uncertainty Propagation")
dx = 100.0
print(f"Initial position uncertainty: ±{dx} m in each direction")
delta_r0 = np.array([dx, dx, dx, 0.0, 0.0, 0.0])
delta_rf = phi @ delta_r0
print(
    f"Final position uncertainty: [{delta_rf[0]:.1f}, {delta_rf[1]:.1f}, {delta_rf[2]:.1f}] m"
)
print(
    f"Uncertainty growth: {np.linalg.norm(delta_rf[0:3]) / np.linalg.norm(delta_r0[0:3]):.1f}x"
)

# Example output:
# Integration with State Transition Matrix
# Initial orbit: 6978.1 km altitude
# Time   Position STM[0,0]  Velocity STM[3,3]  Det(STM)
# ------------------------------------------------------------
#      0s      1.0000        1.0000          1.0000
#    223s      1.0580        1.0564          1.0000
#    594s      1.3993        1.3227          1.0000
#   1007s      2.0473        1.5071          1.0000
#   1346s      2.5985        1.1989          1.0000
#   1530s      2.8009        0.7891          1.0000
#   1849s      2.7780       -0.2849          1.0000
#   2245s      1.6741       -1.8673          1.0000
#   2608s     -0.7226       -2.9191          1.0000
#   2835s     -2.8726       -3.1249          1.0000
#   3091s     -5.7502       -2.8654          1.0000
#   3455s    -10.0598       -1.7595          1.0000
#   3850s    -13.8989       -0.1764          1.0000
#   4169s    -15.4760        0.8634          1.0000
#   4360s    -15.5400        1.2575          1.0000
#   4700s    -13.8708        1.5097          1.0000
#   5114s     -8.9723        1.2931          1.0000
#   5484s     -2.6152        1.0402          1.0000
#   5697s      1.5156        1.0008          1.0000

# Propagation complete! (370 steps)

# Example: Uncertainty Propagation
# Initial position uncertainty: ±100.0 m in each direction
# Final position uncertainty: [357.4, -1684.0, 99.0] m
# Uncertainty growth: 10.0x

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
use brahe::*;
use nalgebra::{DMatrix, DVector};

fn main() {
    // Initialize EOP
    initialize_eop().unwrap();

    // Define two-body dynamics (for integrator - takes optional params)
    let dynamics = |_t: f64, state: &DVector<f64>, _params: Option<&DVector<f64>>| -> DVector<f64> {
        let r = nalgebra::Vector3::new(state[0], state[1], state[2]);
        let v = nalgebra::Vector3::new(state[3], state[4], state[5]);
        let r_norm = r.norm();
        let a = -GM_EARTH / r_norm.powi(3) * r;
        DVector::from_vec(vec![v[0], v[1], v[2], a[0], a[1], a[2]])
    };

    // Same dynamics for Jacobian (without params argument)
    let dynamics_for_jac = |_t: f64, state: &DVector<f64>, _params: Option<&DVector<f64>>| -> DVector<f64> {
        let r = nalgebra::Vector3::new(state[0], state[1], state[2]);
        let v = nalgebra::Vector3::new(state[3], state[4], state[5]);
        let r_norm = r.norm();
        let a = -GM_EARTH / r_norm.powi(3) * r;
        DVector::from_vec(vec![v[0], v[1], v[2], a[0], a[1], a[2]])
    };

    // Create numerical Jacobian for variational equations
    let jacobian = DNumericalJacobian::central(Box::new(dynamics_for_jac))
        .with_adaptive(1.0, 1e-6);

    // Initial orbit (LEO)
    let r0 = DVector::from_vec(vec![R_EARTH + 600e3, 0.0, 0.0]);
    let v0 = DVector::from_vec(vec![0.0, 7.5e3, 0.0]);
    let state0 = DVector::from_vec(vec![r0[0], r0[1], r0[2], v0[0], v0[1], v0[2]]);

    // Initial state transition matrix (identity)
    let phi0 = DMatrix::identity(6, 6);

    println!("Integration with State Transition Matrix");
    println!("Initial orbit: {:.1} km altitude", r0[0] / 1e3);

    // Create integrator with Jacobian
    let config = IntegratorConfig::adaptive(1e-12, 1e-11);
    let integrator = DormandPrince54DIntegrator::with_config(6, Box::new(dynamics), Some(Box::new(jacobian)), None, None, config);

    // Propagate for one orbit period
    let mut t = 0.0;
    let mut state = state0.clone();
    let mut phi = phi0.clone();
    let mut dt: f64 = 60.0;

    // Approximate orbital period
    let period = orbital_period(r0.norm());

    println!("Time   Position STM[0,0]  Velocity STM[3,3]  Det(STM)");
    println!("{}", "-".repeat(60));

    let mut steps = 0;
    while t < period {
        // Propagate state and STM together
        let result = integrator.step_with_varmat(t, state, None, phi, Some(dt.min(period - t)));

        let new_state = result.state;
        let new_phi = result.phi.unwrap();
        let dt_used = result.dt_used;
        let dt_next = result.dt_next;

        t += dt_used;
        state = new_state;
        phi = new_phi;
        dt = dt_next;
        steps += 1;

        // Print progress
        if steps % 20 == 1 {
            let det_phi = phi.clone().lu().determinant();
            println!("{:6.0}s    {:8.4}      {:8.4}        {:8.4}",
                     t, phi[(0, 0)], phi[(3, 3)], det_phi);
        }
    }

    println!("\nPropagation complete! ({} steps)", steps);

    // Example: Map initial position uncertainty to final uncertainty
    println!("\nExample: Uncertainty Propagation");
    let dx = 100.0;
    println!("Initial position uncertainty: ±{} m in each direction", dx);
    let delta_r0 = DVector::from_vec(vec![dx, dx, dx, 0.0, 0.0, 0.0]);
    let delta_rf = &phi * &delta_r0;
    println!("Final position uncertainty: [{:.1}, {:.1}, {:.1}] m",
             delta_rf[0], delta_rf[1], delta_rf[2]);

    let r0_norm = (delta_r0[0].powi(2) + delta_r0[1].powi(2) + delta_r0[2].powi(2)).sqrt();
    let rf_norm = (delta_rf[0].powi(2) + delta_rf[1].powi(2) + delta_rf[2].powi(2)).sqrt();
    println!("Uncertainty growth: {:.1}x", rf_norm / r0_norm);
}

// Example output:
// Integration with State Transition Matrix
// Initial orbit: 6978.1 km altitude
// Time   Position STM[0,0]  Velocity STM[3,3]  Det(STM)
// ------------------------------------------------------------
//      0s      1.0000        1.0000          1.0000
//    223s      1.0580        1.0564          1.0000
//    594s      1.3993        1.3227          1.0000
//   1007s      2.0473        1.5071          1.0000
//   1346s      2.5985        1.1989          1.0000
//   1530s      2.8009        0.7891          1.0000
//   1849s      2.7780       -0.2849          1.0000
//   2245s      1.6741       -1.8673          1.0000
//   2608s     -0.7226       -2.9191          1.0000
//   2835s     -2.8726       -3.1249          1.0000
//   3091s     -5.7502       -2.8654          1.0000
//   3455s    -10.0598       -1.7595          1.0000
//   3850s    -13.8989       -0.1764          1.0000
//   4169s    -15.4760        0.8634          1.0000
//   4360s    -15.5400        1.2575          1.0000
//   4700s    -13.8708        1.5097          1.0000
//   5114s     -8.9723        1.2931          1.0000
//   5484s     -2.6152        1.0402          1.0000
//   5697s      1.5156        1.0008          1.0000

// Propagation complete! (370 steps)

// Example: Uncertainty Propagation
// Initial position uncertainty: ±100.0 m in each direction
// Final position uncertainty: [357.4, -1684.0, 99.0] m
// Uncertainty growth: 10.0x

Equivalence to Direct Perturbation

The power of the STM is that it allows predicting many perturbed trajectories efficiently. Instead of integrating each perturbed initial condition separately, we can integrate the nominal trajectory once (with STM) and map any initial perturbation through the STM.

The following example demonstrates this equivalence:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
import brahe as bh
import numpy as np

# Initialize EOP
bh.initialize_eop()


# Define two-body orbital dynamics
def dynamics(t, state):
    """Two-body point-mass dynamics with Earth gravity."""
    mu = bh.GM_EARTH
    r = state[0:3]
    v = state[3:6]
    r_norm = np.linalg.norm(r)
    a = -mu / r_norm**3 * r
    return np.concatenate([v, a])


# Create numerical Jacobian for variational equations
jacobian = bh.NumericalJacobian.central(dynamics).with_fixed_offset(0.1)

# Configuration for high accuracy
config = bh.IntegratorConfig.adaptive(abs_tol=1e-12, rel_tol=1e-10)

# Create two integrators:
# 1. With Jacobian - propagates STM alongside state
integrator_nominal = bh.RKN1210Integrator(6, dynamics, jacobian=jacobian, config=config)

# 2. Without Jacobian - for direct perturbation integration
integrator_pert = bh.RKN1210Integrator(6, dynamics, config=config)

# Initial state: circular orbit at 500 km altitude
oe0 = np.array([bh.R_EARTH + 500e3, 0.0, 0.0, 0.0, 0.0, 0.0])
state0 = bh.state_koe_to_eci(oe0, bh.AngleFormat.DEGREES)

# Apply 10-meter perturbation in X direction
perturbation = np.array([10.0, 0.0, 0.0, 0.0, 0.0, 0.0])

# Integration parameters
total_time = 100.0  # Total propagation time (seconds)
num_steps = 10  # Number of steps
dt = total_time / num_steps

# Initialize states
state_nominal = state0.copy()
phi = np.eye(6)  # State Transition Matrix starts as identity
state_pert = state0 + perturbation

print("STM vs Direct Perturbation Comparison")
print("=" * 70)
print(f"Initial orbit: {(oe0[0] - bh.R_EARTH) / 1e3:.1f} km altitude (circular)")
print(f"Perturbation: {perturbation[0]:.1f} m in X direction")
print(f"Propagating for {total_time:.0f} seconds in {num_steps} steps\n")
print("Theory: For small δx₀, the perturbed state should satisfy:")
print("        x_pert(t) ≈ x_nominal(t) + Φ(t,t₀)·δx₀\n")
print("Step   Time(s)  ||Error||(m)  Max Component(m)  Relative Error")
print("-" * 70)

t = 0.0
for step in range(num_steps):
    # Propagate nominal trajectory with STM
    new_state_nominal, new_phi, dt_used, _, _ = integrator_nominal.step_with_varmat(
        t, state_nominal, phi, dt
    )

    # Propagate perturbed trajectory directly
    result_pert = integrator_pert.step(t, state_pert, dt)

    # Predict perturbed state using STM: x_pert ≈ x_nominal + Φ·δx₀
    state_pert_predicted = new_state_nominal + new_phi @ perturbation

    # Compute error between STM prediction and direct integration
    error = result_pert.state - state_pert_predicted
    error_norm = np.linalg.norm(error)
    error_max = np.max(np.abs(error))

    # Relative error compared to perturbation magnitude
    relative_error = error_norm / np.linalg.norm(perturbation)

    print(
        f"{step + 1:4d}  {t + dt_used:7.1f}  {error_norm:12.6f}  {error_max:16.6f}  {relative_error:13.6f}"
    )

    # Update for next step
    state_nominal = new_state_nominal
    phi = new_phi
    state_pert = result_pert.state
    t += dt_used

# Example output:
# STM vs Direct Perturbation Comparison
# ======================================================================
# Initial orbit: 500.0 km altitude (circular)
# Perturbation: 10.0 m in X direction
# Propagating for 100 seconds in 10 steps
#
# Theory: For small δx₀, the perturbed state should satisfy:
#         x_pert(t) ≈ x_nominal(t) + Φ(t,t₀)·δx₀
#
# Step   Time(s)  ||Error||(m)  Max Component(m)  Relative Error
# ----------------------------------------------------------------------
#    1     10.0      0.000078          0.000053      0.000008
#    2     20.0      0.000299          0.000211      0.000030
#    3     30.0      0.000627          0.000462      0.000063
#    4     40.0      0.001025          0.000791      0.000103
#    5     50.0      0.001463          0.001176      0.000146
#    6     60.0      0.001919          0.001600      0.000192
#    7     70.0      0.002378          0.002057      0.000238
#    8     80.0      0.002831          0.002539      0.000283
#    9     90.0      0.003271          0.003040      0.000327
#   10    100.0      0.003693          0.003556      0.000369

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
use brahe::*;
use nalgebra::{DMatrix, DVector, SVector};

/// Two-body point-mass dynamics with Earth gravity (for integrator)
fn dynamics(_t: f64, state: &DVector<f64>, _params: Option<&DVector<f64>>) -> DVector<f64> {
    let r = state.rows(0, 3);
    let v = state.rows(3, 3);
    let r_norm = r.norm();
    let a = -GM_EARTH / (r_norm * r_norm * r_norm) * r;

    let mut state_dot = DVector::zeros(6);
    state_dot.rows_mut(0, 3).copy_from(&v);
    state_dot.rows_mut(3, 3).copy_from(&a);
    state_dot
}

/// Two-body dynamics (for Jacobian computation - no params)
fn dynamics_for_jac(_t: f64, state: &DVector<f64>, _params: Option<&DVector<f64>>) -> DVector<f64> {
    let r = state.rows(0, 3);
    let v = state.rows(3, 3);
    let r_norm = r.norm();
    let a = -GM_EARTH / (r_norm * r_norm * r_norm) * r;

    let mut state_dot = DVector::zeros(6);
    state_dot.rows_mut(0, 3).copy_from(&v);
    state_dot.rows_mut(3, 3).copy_from(&a);
    state_dot
}

fn main() {
    // Initialize EOP
    initialize_eop().unwrap();

    // Create numerical Jacobian for variational equations
    let jacobian = DNumericalJacobian::central(Box::new(dynamics_for_jac))
        .with_fixed_offset(0.1);

    // Configuration for high accuracy
    let config = IntegratorConfig::adaptive(1e-12, 1e-10);

    // Create two integrators:
    // 1. With Jacobian - propagates STM alongside state
    let integrator_nominal = RKN1210DIntegrator::with_config(
        6,
        Box::new(dynamics),
        Some(Box::new(jacobian)),
        None,
        None,
        config.clone(),
    );

    // 2. Without Jacobian - for direct perturbation integration
    let integrator_pert = RKN1210DIntegrator::with_config(
        6,
        Box::new(dynamics),
        None,
        None,
        None,
        config,
    );

    // Initial state: circular orbit at 500 km altitude
    let oe0 = SVector::<f64, 6>::new(R_EARTH + 500e3, 0.0, 0.0, 0.0, 0.0, 0.0);
    let state0_static = state_koe_to_eci(oe0, AngleFormat::Degrees);
    let state0 = DVector::from_vec(state0_static.as_slice().to_vec());

    // Apply 10-meter perturbation in X direction
    let perturbation = DVector::from_vec(vec![10.0, 0.0, 0.0, 0.0, 0.0, 0.0]);

    // Integration parameters
    let total_time = 100.0; // Total propagation time (seconds)
    let num_steps = 10;     // Number of steps
    let dt = total_time / num_steps as f64;

    // Initialize states
    let mut state_nominal = state0.clone();
    let mut phi = DMatrix::identity(6, 6); // State Transition Matrix starts as identity
    let mut state_pert = &state0 + &perturbation;
    let mut t = 0.0;

    println!("STM vs Direct Perturbation Comparison");
    println!("{}", "=".repeat(70));
    println!("Initial orbit: {:.1} km altitude (circular)", (oe0[0] - R_EARTH) / 1e3);
    println!("Perturbation: {:.1} m in X direction", perturbation[0]);
    println!("Propagating for {:.0} seconds in {} steps\n", total_time, num_steps);
    println!("Theory: For small δx₀, the perturbed state should satisfy:");
    println!("        x_pert(t) ≈ x_nominal(t) + Φ(t,t₀)·δx₀\n");
    println!("Step   Time(s)  ||Error||(m)  Max Component(m)  Relative Error");
    println!("{}", "-".repeat(70));

    for step in 0..num_steps {
        // Propagate nominal trajectory with STM
        let result_nominal = integrator_nominal.step_with_varmat(
            t,
            state_nominal.clone(),
            None,
            phi.clone(),
            Some(dt),
        );
        let new_state_nominal = result_nominal.state;
        let new_phi = result_nominal.phi.unwrap();
        let dt_used = result_nominal.dt_used;

        // Propagate perturbed trajectory directly
        let result_pert = integrator_pert.step(
            t,
            state_pert.clone(),
            None,
            Some(dt),
        );

        // Predict perturbed state using STM: x_pert ≈ x_nominal + Φ·δx₀
        let state_pert_predicted = &new_state_nominal + &new_phi * &perturbation;

        // Compute error between STM prediction and direct integration
        let error = &result_pert.state - &state_pert_predicted;
        let error_norm = error.norm();
        let error_max = error.abs().max();

        // Relative error compared to perturbation magnitude
        let relative_error = error_norm / perturbation.norm();

        println!(
            "{:4}  {:7.1}  {:12.6}  {:16.6}  {:13.6}",
            step + 1,
            t + dt_used,
            error_norm,
            error_max,
            relative_error
        );

        // Update for next step
        state_nominal = new_state_nominal;
        phi = new_phi;
        state_pert = result_pert.state;
        t += dt_used;
    }
}

// Example output:
// STM vs Direct Perturbation Comparison
// ======================================================================
// Initial orbit: 500.0 km altitude (circular)
// Perturbation: 10.0 m in X direction
// Propagating for 100 seconds in 10 steps
//
// Theory: For small δx₀, the perturbed state should satisfy:
//         x_pert(t) ≈ x_nominal(t) + Φ(t,t₀)·δx₀
//
// Step   Time(s)  ||Error||(m)  Max Component(m)  Relative Error
// ----------------------------------------------------------------------
//    1     10.0      0.000078          0.000053      0.000008
//    2     20.0      0.000299          0.000211      0.000030
//    3     30.0      0.000627          0.000462      0.000063
//    4     40.0      0.001025          0.000791      0.000103
//    5     50.0      0.001463          0.001176      0.000146
//    6     60.0      0.001919          0.001600      0.000192
//    7     70.0      0.002378          0.002057      0.000238
//    8     80.0      0.002831          0.002539      0.000283
//    9     90.0      0.003271          0.003040      0.000327
//   10    100.0      0.003693          0.003556      0.000369

See Also