Skip to content

Extending Spacecraft State

The NumericalOrbitPropagator supports extending state vectors beyond the standard 6-element orbital state, enabling modeling of additional state variables and dynamics like propellant mass, battery charge, or attiude alongside orbital dynamics. This is achieved through the additional_dynamics function.

State Extension Approach

To extend the state vector with NumericalOrbitPropagator:

  1. Define an extended state vector (e.g., 7 elements: [pos, vel, mass])
  2. Implement an additional_dynamics function that returns a full state-sized derivative vector, where the first 6 elements are zeros (orbital dynamics handled by the force model) and the remaining elements contain derivatives for the extended state
  3. Optionally provide a control_input function for thrust accelerations
  4. Create the propagator with these functions

The key advantage of using NumericalOrbitPropagator is that orbital dynamics (gravity, drag, SRP, etc.) are handled automatically by the force model configuration, while your additional_dynamics function adds derivatives for the extended state elements.

To showcase how to extend the spacecraft state, we present an example of tracking propellant mass during a thrust maneuver below.

Mass Tracking Example

One common extension is tracking propellant mass during the mission. To model propelant mass we augment the state vector from 6 to 7 elements, by adding mass \(m\) as the 7th element:

\[\mathbf{x} = [x, y, z, v_x, v_y, v_z, m]^T\]

Mass Flow Dynamics

We model mass flow rate during thrust as:

\[\dot{m} = -\frac{F}{I_{sp} \cdot g_0}\]

where:

  • \(F\) is thrust force (N)
  • \(I_{sp}\) is specific impulse (s)
  • \(g_0\) is standard gravity (9.80665 m/s²)

Implementation with NumericalOrbitPropagator

Both additional_dynamics and control_input functions return full state-sized vectors. The propagator adds these to the orbital dynamics computed from the force model.

The additional_dynamics function returns a state-sized vector with derivatives for extended elements:

1
2
3
4
5
6
def additional_dynamics(t, state, params):
    """Return full state-sized vector with mass rate."""
    dx = np.zeros(len(state))  # Full state size
    if burn_start <= t < burn_end:
        dx[6] = -mass_flow_rate  # dm/dt = -F/(Isp*g0)
    return dx

1
2
3
4
5
6
7
8
9
let additional_dynamics: bh::DAdditionalDynamics = Some(Box::new(
    |t: f64, state: &na::DVector<f64>, _params: Option<&na::DVector<f64>>| {
        let mut dx = na::DVector::zeros(state.len());
        if burn_start <= t && t < burn_end {
            dx[6] = -mass_flow_rate;  // dm/dt = -F/(Isp*g0)
        }
        dx
    },
));

The control_input function returns a state-sized vector with acceleration in indices 3-5:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
def control_input(t, state, params):
    """Return full state-sized vector with thrust acceleration."""
    dx = np.zeros(len(state))
    if burn_start <= t < burn_end:
        mass = state[6]  # Access mass from extended state
        vel = state[3:6]
        v_hat = vel / np.linalg.norm(vel)  # Prograde direction
        acc = (thrust_force / mass) * v_hat
        dx[3:6] = acc  # Add to velocity derivatives
    return dx

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
let control_input: bh::DControlInput = Some(Box::new(
    |t: f64, state: &na::DVector<f64>, _params: Option<&na::DVector<f64>>| {
        let mut dx = na::DVector::zeros(state.len());
        if burn_start <= t && t < burn_end {
            let mass = state[6];  // Access mass from extended state
            let vel = na::Vector3::new(state[3], state[4], state[5]);
            let v_hat = vel.normalize();  // Prograde direction
            let acc = (thrust_force / mass) * v_hat;
            dx[3] = acc[0];
            dx[4] = acc[1];
            dx[5] = acc[2];
        }
        dx
    },
));

Complete Example

  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
import numpy as np
import brahe as bh

# Initialize EOP data
bh.initialize_eop()

# Create initial epoch
epoch = bh.Epoch.from_datetime(2024, 1, 1, 12, 0, 0.0, 0.0, bh.TimeSystem.UTC)

# Initial orbital elements and state
oe = np.array([bh.R_EARTH + 500e3, 0.01, 45.0, 15.0, 30.0, 45.0])
orbital_state = bh.state_koe_to_eci(oe, bh.AngleFormat.DEGREES)

# Extended state: [x, y, z, vx, vy, vz, mass]
initial_mass = 1000.0  # kg
initial_state = np.concatenate([orbital_state, [initial_mass]])

# Thruster parameters
thrust_force = 10.0  # N
specific_impulse = 300.0  # s
g0 = 9.80665  # m/s^2
mass_flow_rate = thrust_force / (specific_impulse * g0)  # kg/s

# Timing parameters
pre_burn_coast = 300.0  # 5 minutes coast before burn
burn_duration = 600.0  # 10 minutes burn
post_burn_coast = 600.0  # 10 minutes coast after burn
burn_start = pre_burn_coast
burn_end = pre_burn_coast + burn_duration
total_time = pre_burn_coast + burn_duration + post_burn_coast

# Spacecraft parameters for force model [mass, drag_area, Cd, srp_area, Cr]
params = np.array([initial_mass, 2.0, 2.2, 2.0, 1.3])

print("Thruster parameters:")
print(f"  Thrust: {thrust_force} N")
print(f"  Isp: {specific_impulse} s")
print(f"  Mass flow rate: {mass_flow_rate * 1000:.2f} g/s")
print(f"  Burn duration: {burn_duration} s")
print(f"  Burn window: {burn_start} - {burn_end} s")
print(f"  Expected fuel consumption: {mass_flow_rate * burn_duration:.2f} kg")


# Define additional dynamics for mass tracking
def additional_dynamics(t, state, params):
    """
    Return full state derivative vector with contributions for extended state.
    State: [x, y, z, vx, vy, vz, mass] - return same-sized vector.
    Elements 0-5 should be zero (orbital dynamics handled by force model).
    """
    dx = np.zeros(len(state))
    if burn_start <= t < burn_end:
        dx[6] = -mass_flow_rate  # dm/dt = -F/(Isp*g0)
    return dx


# Define control input for thrust acceleration
def control_input(t, state, params):
    """
    Return full state derivative with acceleration contributions.
    Returns state-sized vector with acceleration in indices 3-5.
    """
    dx = np.zeros(len(state))
    if burn_start <= t < burn_end:
        mass = state[6]  # Access mass from extended state
        vel = state[3:6]
        v_hat = vel / np.linalg.norm(vel)  # Prograde direction
        acc = (thrust_force / mass) * v_hat  # Thrust acceleration
        dx[3:6] = acc  # Add to velocity derivatives
    return dx


# Create propagator with two-body dynamics (no drag/SRP for clean mass tracking)
force_config = bh.ForceModelConfig.two_body()
prop_config = bh.NumericalPropagationConfig.default()

prop = bh.NumericalOrbitPropagator(
    epoch,
    initial_state,
    prop_config,
    force_config,
    params=params,
    additional_dynamics=additional_dynamics,
    control_input=control_input,
)

print("\nInitial state:")
print(f"  Mass: {initial_mass:.1f} kg")
print(f"  Semi-major axis: {oe[0] / 1e3:.1f} km")

# Propagate through pre-burn coast, burn, and post-burn coast
prop.propagate_to(epoch + total_time)

# Check final state
final_state = prop.current_state()
final_mass = final_state[6]
fuel_consumed = initial_mass - final_mass

# Compute final orbital elements
final_orbital_state = final_state[:6]
final_koe = bh.state_eci_to_koe(final_orbital_state, bh.AngleFormat.DEGREES)

print("\nFinal state:")
print(f"  Mass: {final_mass:.1f} kg")
print(f"  Fuel consumed: {fuel_consumed:.2f} kg")
print(f"  Semi-major axis: {final_koe[0] / 1e3:.1f} km")
print(f"  Delta-a: {(final_koe[0] - oe[0]) / 1e3:.1f} km")

# Verify Tsiolkovsky equation
delta_v_expected = specific_impulse * g0 * np.log(initial_mass / final_mass)
print("\nTsiolkovsky verification:")
print(f"  Expected delta-v: {delta_v_expected:.1f} m/s")

# Validate
expected_fuel = mass_flow_rate * burn_duration
assert abs(fuel_consumed - expected_fuel) < 0.1  # Within 0.1 kg
assert final_mass < initial_mass
assert final_koe[0] > oe[0]  # Orbit raised

print("\nExample validated successfully!")

  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
use brahe as bh;
use bh::integrators::traits::DStateDynamics;
use bh::propagators::{DNumericalOrbitPropagator, ForceModelConfig, NumericalPropagationConfig};
use bh::traits::DStatePropagator;
use nalgebra as na;

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

    // Create initial epoch
    let epoch = bh::Epoch::from_datetime(2024, 1, 1, 12, 0, 0.0, 0.0, bh::TimeSystem::UTC);

    // Initial orbital elements and state
    let oe = na::SVector::<f64, 6>::new(bh::R_EARTH + 500e3, 0.01, 45.0, 0.0, 0.0, 0.0);
    let orbital_state = bh::state_koe_to_eci(oe, bh::AngleFormat::Degrees);

    // Extended state: [x, y, z, vx, vy, vz, mass]
    let initial_mass = 1000.0; // kg
    let mut initial_state = na::DVector::zeros(7);
    for i in 0..6 {
        initial_state[i] = orbital_state[i];
    }
    initial_state[6] = initial_mass;

    // Thruster parameters
    let thrust_force = 10.0; // N
    let specific_impulse = 300.0; // s
    let g0 = 9.80665; // m/s^2
    let mass_flow_rate = thrust_force / (specific_impulse * g0); // kg/s

    // Timing parameters
    let pre_burn_coast = 300.0; // 5 minutes coast before burn
    let burn_duration = 600.0; // 10 minutes burn
    let post_burn_coast = 600.0; // 10 minutes coast after burn
    let burn_start = pre_burn_coast;
    let burn_end = pre_burn_coast + burn_duration;
    let total_time = pre_burn_coast + burn_duration + post_burn_coast;

    println!("Thruster parameters:");
    println!("  Thrust: {} N", thrust_force);
    println!("  Isp: {} s", specific_impulse);
    println!("  Mass flow rate: {:.2} g/s", mass_flow_rate * 1000.0);
    println!("  Burn duration: {} s", burn_duration);
    println!("  Burn window: {} - {} s", burn_start, burn_end);
    println!(
        "  Expected fuel consumption: {:.2} kg",
        mass_flow_rate * burn_duration
    );

    // Additional dynamics for mass tracking
    // Returns full state-sized vector with mass derivative
    let additional_dynamics: DStateDynamics = Box::new(move |t, state, _params| {
        let mut dx = na::DVector::zeros(state.len());
        if burn_start <= t && t < burn_end {
            dx[6] = -mass_flow_rate; // dm/dt = -F/(Isp*g0)
        }
        dx
    });

    // Control input for thrust acceleration
    // Returns full state-sized vector with acceleration in indices 3-5
    let control_input = Some(Box::new(move |t: f64, state: &na::DVector<f64>, _params: Option<&na::DVector<f64>>| {
        let mut dx = na::DVector::zeros(state.len());
        if burn_start <= t && t < burn_end {
            let mass = state[6];
            let vel = state.fixed_rows::<3>(3);
            let v_hat = vel.normalize();
            let acc = (thrust_force / mass) * v_hat;
            dx[3] = acc[0];
            dx[4] = acc[1];
            dx[5] = acc[2];
        }
        dx
    }) as Box<dyn Fn(f64, &na::DVector<f64>, Option<&na::DVector<f64>>) -> na::DVector<f64> + Send + Sync>);

    // Create propagator with two-body dynamics (no drag/SRP for clean mass tracking)
    let mut prop = DNumericalOrbitPropagator::new(
        epoch,
        initial_state.clone(),
        NumericalPropagationConfig::default(),
        ForceModelConfig::earth_gravity(),
        None, // params
        Some(additional_dynamics),
        control_input,
        None, // initial_covariance
    )
    .unwrap();

    println!("\nInitial state:");
    println!("  Mass: {:.1} kg", initial_mass);
    println!("  Semi-major axis: {:.1} km", oe[0] / 1e3);

    // Propagate through pre-burn coast, burn, and post-burn coast
    prop.propagate_to(epoch + total_time);

    // Check final state
    let final_state = prop.current_state();
    let final_mass = final_state[6];
    let fuel_consumed = initial_mass - final_mass;

    // Compute final orbital elements
    let final_orbital_state =
        na::SVector::<f64, 6>::from_column_slice(&final_state.as_slice()[..6]);
    let final_koe = bh::state_eci_to_koe(final_orbital_state, bh::AngleFormat::Degrees);

    println!("\nFinal state:");
    println!("  Mass: {:.1} kg", final_mass);
    println!("  Fuel consumed: {:.2} kg", fuel_consumed);
    println!("  Semi-major axis: {:.1} km", final_koe[0] / 1e3);
    println!("  Delta-a: {:.1} km", (final_koe[0] - oe[0]) / 1e3);

    // Verify Tsiolkovsky equation
    let delta_v_expected = specific_impulse * g0 * (initial_mass / final_mass).ln();
    println!("\nTsiolkovsky verification:");
    println!("  Expected delta-v: {:.1} m/s", delta_v_expected);

    // Validate
    let expected_fuel = mass_flow_rate * burn_duration;
    assert!((fuel_consumed - expected_fuel).abs() < 0.1); // Within 0.1 kg
    assert!(final_mass < initial_mass);
    assert!(final_koe[0] > oe[0]); // Orbit raised

    println!("\nExample validated successfully!");
}

Orbital Elements Evolution

The following plot shows how orbital elements evolve during the thrust maneuver:

Plot Source
mass_tracking_elements.py
  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
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

# Add plots directory to path for importing brahe_theme
sys.path.insert(0, str(pathlib.Path(__file__).parent.parent.parent))
from brahe_theme import save_themed_html, get_theme_colors

# Configuration
SCRIPT_NAME = pathlib.Path(__file__).stem
OUTDIR = pathlib.Path(os.getenv("BRAHE_FIGURE_OUTPUT_DIR", "./docs/figures/"))
os.makedirs(OUTDIR, exist_ok=True)

# Initialize EOP data
bh.initialize_eop()

# Create initial epoch and state
epoch = bh.Epoch.from_datetime(2024, 1, 1, 12, 0, 0.0, 0.0, bh.TimeSystem.UTC)
oe = np.array([bh.R_EARTH + 500e3, 0.01, 45.0, 15.0, 30.0, 45.0])
orbital_state = bh.state_koe_to_eci(oe, bh.AngleFormat.DEGREES)

# Extended state: [x, y, z, vx, vy, vz, mass]
initial_mass = 1000.0  # kg
initial_state = np.concatenate([orbital_state, [initial_mass]])

# Thruster parameters
thrust_force = 10.0  # N
specific_impulse = 300.0  # s
g0 = 9.80665  # m/s^2
mass_flow_rate = thrust_force / (specific_impulse * g0)  # kg/s

# Timing parameters
pre_burn_coast = 300.0  # 5 minutes coast before burn
burn_duration = 600.0  # 10 minutes burn
post_burn_coast = 600.0  # 10 minutes coast after burn
burn_start = pre_burn_coast
burn_end = pre_burn_coast + burn_duration
total_time = pre_burn_coast + burn_duration + post_burn_coast

# Spacecraft parameters for force model
params = np.array([initial_mass, 2.0, 2.2, 2.0, 1.3])


# Define additional dynamics for mass tracking
def additional_dynamics(t, state, params):
    dx = np.zeros(len(state))
    if burn_start <= t < burn_end:
        dx[6] = -mass_flow_rate
    return dx


# Define control input for thrust acceleration
def control_input(t, state, params):
    dx = np.zeros(len(state))
    if burn_start <= t < burn_end:
        mass = state[6]
        vel = state[3:6]
        v_hat = vel / np.linalg.norm(vel)
        acc = (thrust_force / mass) * v_hat
        dx[3:6] = acc
    return dx


# Create propagator with two-body dynamics
force_config = bh.ForceModelConfig.two_body()
prop_config = bh.NumericalPropagationConfig.default()

prop = bh.NumericalOrbitPropagator(
    epoch,
    initial_state,
    prop_config,
    force_config,
    params=params,
    additional_dynamics=additional_dynamics,
    control_input=control_input,
)

# Propagate and collect orbital elements over time
prop.propagate_to(epoch + total_time)

# Sample trajectory using trajectory interpolation
traj = prop.trajectory
times = []
a_vals = []  # Semi-major axis (km)
e_vals = []  # Eccentricity
i_vals = []  # Inclination (deg)
raan_vals = []  # RAAN (deg)
argp_vals = []  # Argument of periapsis (deg)
ma_vals = []  # Mean anomaly (deg)

dt = 10.0  # 10 second samples
t = 0.0
while t <= total_time:
    current_epoch = epoch + t
    try:
        state = traj.interpolate(current_epoch)
        # Convert to orbital elements
        orbital_state_6d = state[:6]
        koe = bh.state_eci_to_koe(orbital_state_6d, bh.AngleFormat.DEGREES)

        times.append(t / 60.0)  # Convert to minutes
        a_vals.append((koe[0] - bh.R_EARTH) / 1e3)  # Altitude in km
        e_vals.append(koe[1])
        i_vals.append(koe[2])
        raan_vals.append(koe[3])
        argp_vals.append(koe[4])
        ma_vals.append(koe[5])  # Mean anomaly (koe[5] is mean anomaly)
    except RuntimeError:
        pass  # Skip if interpolation fails

    t += dt


def create_figure(theme):
    colors = get_theme_colors(theme)

    # Create 2x3 subplots
    fig = make_subplots(
        rows=2,
        cols=3,
        subplot_titles=(
            "Semi-major Axis (Altitude)",
            "Eccentricity",
            "Inclination",
            "RAAN",
            "Arg. Periapsis",
            "Mean Anomaly",
        ),
        vertical_spacing=0.15,
        horizontal_spacing=0.08,
    )

    # Semi-major axis (altitude)
    fig.add_trace(
        go.Scatter(
            x=times, y=a_vals, mode="lines", line=dict(color=colors["primary"], width=2)
        ),
        row=1,
        col=1,
    )
    fig.update_yaxes(title_text="Altitude (km)", row=1, col=1)

    # Eccentricity
    fig.add_trace(
        go.Scatter(
            x=times,
            y=e_vals,
            mode="lines",
            line=dict(color=colors["secondary"], width=2),
        ),
        row=1,
        col=2,
    )
    fig.update_yaxes(title_text="e", range=[0, 0.1], row=1, col=2)

    # Inclination
    fig.add_trace(
        go.Scatter(
            x=times, y=i_vals, mode="lines", line=dict(color=colors["accent"], width=2)
        ),
        row=1,
        col=3,
    )
    fig.update_yaxes(title_text="i (deg)", range=[0, 90], row=1, col=3)

    # RAAN
    fig.add_trace(
        go.Scatter(
            x=times,
            y=raan_vals,
            mode="lines",
            line=dict(color=colors["primary"], width=2),
        ),
        row=2,
        col=1,
    )
    fig.update_yaxes(title_text="RAAN (deg)", range=[0, 360], row=2, col=1)

    # Argument of periapsis
    fig.add_trace(
        go.Scatter(
            x=times,
            y=argp_vals,
            mode="lines",
            line=dict(color=colors["secondary"], width=2),
        ),
        row=2,
        col=2,
    )
    fig.update_yaxes(title_text="\u03c9 (deg)", range=[0, 360], row=2, col=2)

    # Mean anomaly
    fig.add_trace(
        go.Scatter(
            x=times,
            y=ma_vals,
            mode="lines",
            line=dict(color=colors["accent"], width=2),
        ),
        row=2,
        col=3,
    )
    fig.update_yaxes(title_text="M (deg)", range=[0, 360], row=2, col=3)

    # Add burn start and end indicators to all subplots
    burn_start_min = burn_start / 60.0
    burn_end_min = burn_end / 60.0
    for row in [1, 2]:
        for col in [1, 2, 3]:
            # Burn start indicator
            fig.add_vline(
                x=burn_start_min,
                line_dash="dot",
                line_color=colors["accent"],
                line_width=1,
                row=row,
                col=col,
            )
            # Burn end indicator
            fig.add_vline(
                x=burn_end_min,
                line_dash="dot",
                line_color=colors["error"],
                line_width=1,
                row=row,
                col=col,
            )

    # Update x-axis labels for bottom row
    for col in [1, 2, 3]:
        fig.update_xaxes(title_text="Time (min)", row=2, col=col)

    fig.update_layout(
        title="Orbital Elements During Prograde Thrust (10 N, 10 min burn)",
        showlegend=False,
        height=500,
        margin=dict(l=60, r=40, t=80, b=60),
    )

    return fig


# Save themed HTML files
light_path, dark_path = save_themed_html(create_figure, OUTDIR / SCRIPT_NAME)
print(f"Generated {light_path}")
print(f"Generated {dark_path}")

Mass Depletion Profile

The mass decreases linearly during the thrust phase:

Plot Source
mass_tracking_mass.py
  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
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
import numpy as np
import plotly.graph_objects as go

# Add plots directory to path for importing brahe_theme
sys.path.insert(0, str(pathlib.Path(__file__).parent.parent.parent))
from brahe_theme import save_themed_html, get_theme_colors

# Configuration
SCRIPT_NAME = pathlib.Path(__file__).stem
OUTDIR = pathlib.Path(os.getenv("BRAHE_FIGURE_OUTPUT_DIR", "./docs/figures/"))
os.makedirs(OUTDIR, exist_ok=True)

# Initialize EOP data
bh.initialize_eop()

# Create initial epoch and state
epoch = bh.Epoch.from_datetime(2024, 1, 1, 12, 0, 0.0, 0.0, bh.TimeSystem.UTC)
oe = np.array([bh.R_EARTH + 500e3, 0.01, 45.0, 15.0, 30.0, 45.0])
orbital_state = bh.state_koe_to_eci(oe, bh.AngleFormat.DEGREES)

# Extended state: [x, y, z, vx, vy, vz, mass]
initial_mass = 1000.0  # kg
initial_state = np.concatenate([orbital_state, [initial_mass]])

# Thruster parameters
thrust_force = 10.0  # N
specific_impulse = 300.0  # s
g0 = 9.80665  # m/s^2
mass_flow_rate = thrust_force / (specific_impulse * g0)  # kg/s

# Timing parameters
pre_burn_coast = 300.0  # 5 minutes coast before burn
burn_duration = 600.0  # 10 minutes burn
post_burn_coast = 600.0  # 10 minutes coast after burn
burn_start = pre_burn_coast
burn_end = pre_burn_coast + burn_duration
total_time = pre_burn_coast + burn_duration + post_burn_coast

# Spacecraft parameters for force model
params = np.array([initial_mass, 2.0, 2.2, 2.0, 1.3])


# Define additional dynamics for mass tracking
def additional_dynamics(t, state, params):
    dx = np.zeros(len(state))
    if burn_start <= t < burn_end:
        dx[6] = -mass_flow_rate
    return dx


# Define control input for thrust acceleration
def control_input(t, state, params):
    dx = np.zeros(len(state))
    if burn_start <= t < burn_end:
        mass = state[6]
        vel = state[3:6]
        v_hat = vel / np.linalg.norm(vel)
        acc = (thrust_force / mass) * v_hat
        dx[3:6] = acc
    return dx


# Create propagator with two-body dynamics
force_config = bh.ForceModelConfig.two_body()
prop_config = bh.NumericalPropagationConfig.default()

prop = bh.NumericalOrbitPropagator(
    epoch,
    initial_state,
    prop_config,
    force_config,
    params=params,
    additional_dynamics=additional_dynamics,
    control_input=control_input,
)

# Propagate and collect mass over time
prop.propagate_to(epoch + total_time)

# Sample trajectory using trajectory interpolation
traj = prop.trajectory
times = []
mass_vals = []
thrust_active = []

dt = 5.0  # 5 second samples
t = 0.0
while t <= total_time:
    current_epoch = epoch + t
    try:
        state = traj.interpolate(current_epoch)
        times.append(t / 60.0)  # Convert to minutes
        mass_vals.append(state[6])
        thrust_active.append(1 if burn_start <= t < burn_end else 0)
    except RuntimeError:
        pass  # Skip if interpolation fails

    t += dt

# Compute expected values
expected_fuel = mass_flow_rate * burn_duration
final_mass_expected = initial_mass - expected_fuel
delta_v_expected = specific_impulse * g0 * np.log(initial_mass / final_mass_expected)


def create_figure(theme):
    colors = get_theme_colors(theme)

    fig = go.Figure()

    # Mass profile
    fig.add_trace(
        go.Scatter(
            x=times,
            y=mass_vals,
            mode="lines",
            name="Spacecraft Mass",
            line=dict(color=colors["primary"], width=3),
        )
    )

    # Initial mass reference
    fig.add_hline(
        y=initial_mass,
        line_dash="dot",
        line_color="gray",
        annotation_text=f"Initial: {initial_mass:.0f} kg",
        annotation_position="top right",
    )

    # Final mass reference
    fig.add_hline(
        y=final_mass_expected,
        line_dash="dot",
        line_color="gray",
        annotation_text=f"Final: {final_mass_expected:.1f} kg",
        annotation_position="bottom right",
    )

    # Thrust phase shading
    burn_start_min = burn_start / 60.0
    burn_end_min = burn_end / 60.0
    fig.add_vrect(
        x0=burn_start_min,
        x1=burn_end_min,
        fillcolor=colors["secondary"],
        opacity=0.1,
        layer="below",
        line_width=0,
        annotation_text="Thrust On",
        annotation_position="top left",
    )

    # Burn start indicator
    fig.add_vline(
        x=burn_start_min,
        line_dash="dash",
        line_color=colors["accent"],
        line_width=2,
        annotation_text="Burn Start",
        annotation_position="top left",
    )

    # Burn end indicator
    fig.add_vline(
        x=burn_end_min,
        line_dash="dash",
        line_color=colors["error"],
        line_width=2,
        annotation_text="Burn End",
        annotation_position="top right",
    )

    fig.update_layout(
        title=f"Mass Depletion Profile (F={thrust_force} N, Isp={specific_impulse} s)",
        xaxis_title="Time (min)",
        yaxis_title="Mass (kg)",
        showlegend=False,
        height=500,
        margin=dict(l=60, r=40, t=80, b=60),
    )

    # Add annotation with summary
    summary_text = (
        f"Fuel consumed: {expected_fuel:.2f} kg<br>"
        f"Mass flow rate: {mass_flow_rate * 1000:.2f} g/s<br>"
        f"Expected \u0394v: {delta_v_expected:.1f} m/s"
    )
    fig.add_annotation(
        x=0.98,
        y=0.5,
        xref="paper",
        yref="paper",
        text=summary_text,
        showarrow=False,
        font=dict(size=11),
        align="right",
        bordercolor="gray",
        borderwidth=1,
        borderpad=4,
        bgcolor="white" if theme == "light" else "#1e1e1e",
        opacity=0.9,
    )

    return fig


# Save themed HTML files
light_path, dark_path = save_themed_html(create_figure, OUTDIR / SCRIPT_NAME)
print(f"Generated {light_path}")
print(f"Generated {dark_path}")

Tsiolkovsky Verification

The mass ratio determines achievable \(\Delta v\):

\[\Delta v = I_{sp} \cdot g_0 \cdot \ln\left(\frac{m_0}{m_f}\right)\]

This provides a useful validation check for mass tracking implementations.

Battery Tracking Example

Another common extension is tracking battery state of charge during eclipse and sunlit periods. This models solar panel charging using the conical shadow model for accurate illumination calculation.

We augment the state vector with battery energy \(E_{bat}\) in Watt-hours:

\[\mathbf{x} = [x, y, z, v_x, v_y, v_z, E_{bat}]^T\]

Power Balance Dynamics

The battery state of charge changes based on the power balance:

\[\dot{E}_{bat} = \nu \cdot P_{solar} - P_{load}\]

where:

  • \(\nu\) is the illumination fraction (0 = full shadow, 1 = full sunlight)
  • \(P_{solar}\) is the solar panel output when fully illuminated (W)
  • \(P_{load}\) is the spacecraft power consumption (W)

Implementation

The additional_dynamics function computes the illumination at each timestep using eclipse_conical:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
def additional_dynamics(t, state, params):
    """Battery dynamics with eclipse-aware solar charging."""
    dx = np.zeros(len(state))
    r_eci = state[:3]

    # Get sun position at current epoch
    current_epoch = epoch + t
    r_sun = bh.sun_position(current_epoch)

    # Get illumination fraction (0 = umbra, 0-1 = penumbra, 1 = sunlit)
    illumination = bh.eclipse_conical(r_eci, r_sun)

    # Battery dynamics (Wh/s = W / 3600)
    power_in = illumination * solar_panel_power  # W
    power_out = load_power  # W
    dx[6] = (power_in - power_out) / 3600.0  # Wh/s

    return dx

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
let additional_dynamics: DStateDynamics = Box::new(move |t, state, _params| {
    let mut dx = na::DVector::zeros(state.len());
    let r_eci = na::Vector3::new(state[0], state[1], state[2]);

    // Get sun position at current epoch
    let current_epoch = *epoch_ref + t;
    let r_sun = bh::sun_position(current_epoch);

    // Get illumination fraction (0 = umbra, 0-1 = penumbra, 1 = sunlit)
    let illumination = bh::eclipse_conical(r_eci, r_sun);

    // Battery dynamics (Wh/s = W / 3600)
    let power_in = illumination * solar_panel_power;  // W
    let power_out = load_power;  // W
    dx[6] = (power_in - power_out) / 3600.0;  // Wh/s

    dx
});

Complete Example

  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
import numpy as np
import brahe as bh

# Initialize EOP data
bh.initialize_eop()

# Create initial epoch
epoch = bh.Epoch.from_datetime(2024, 6, 21, 12, 0, 0.0, 0.0, bh.TimeSystem.UTC)

# Initial orbital elements and state - LEO orbit
oe = np.array([bh.R_EARTH + 500e3, 0.01, 45.0, 15.0, 30.0, 45.0])
orbital_state = bh.state_koe_to_eci(oe, bh.AngleFormat.DEGREES)

# Extended state: [x, y, z, vx, vy, vz, battery_charge]
battery_capacity = 100.0  # Wh
initial_charge = 80.0  # Wh (80% SOC)
initial_state = np.concatenate([orbital_state, [initial_charge]])

# Power system parameters
solar_panel_power = 50.0  # W (when fully illuminated)
load_power = 30.0  # W (continuous consumption)

print("Power system parameters:")
print(f"  Battery capacity: {battery_capacity} Wh")
print(
    f"  Initial charge: {initial_charge} Wh ({100 * initial_charge / battery_capacity:.0f}% SOC)"
)
print(f"  Solar panel power: {solar_panel_power} W")
print(f"  Load power: {load_power} W")
print(f"  Net charging rate (sunlit): {solar_panel_power - load_power} W")
print(f"  Net discharge rate (eclipse): {load_power} W")

# Spacecraft parameters for force model [mass, drag_area, Cd, srp_area, Cr]
params = np.array([500.0, 2.0, 2.2, 2.0, 1.3])


# Define additional dynamics for battery tracking
def additional_dynamics(t, state, params):
    """
    Return full state derivative vector with battery charge dynamics.
    State: [x, y, z, vx, vy, vz, battery_charge] - return same-sized vector.
    Elements 0-5 should be zero (orbital dynamics handled by force model).
    """
    dx = np.zeros(len(state))
    r_eci = state[:3]

    # Get sun position at current epoch
    current_epoch = epoch + t
    r_sun = bh.sun_position(current_epoch)

    # Get illumination fraction (0 = umbra, 0-1 = penumbra, 1 = sunlit)
    illumination = bh.eclipse_conical(r_eci, r_sun)

    # Battery dynamics (Wh/s = W / 3600)
    power_in = illumination * solar_panel_power  # W
    power_out = load_power  # W
    charge_rate = (power_in - power_out) / 3600.0  # Wh/s

    # Apply battery limits (0 to capacity)
    charge = state[6]
    if charge >= battery_capacity and charge_rate > 0:
        charge_rate = 0.0  # Battery full
    elif charge <= 0 and charge_rate < 0:
        charge_rate = 0.0  # Battery empty

    dx[6] = charge_rate
    return dx


# Create propagator with two-body dynamics
force_config = bh.ForceModelConfig.two_body()
prop_config = bh.NumericalPropagationConfig.default()

prop = bh.NumericalOrbitPropagator(
    epoch,
    initial_state,
    prop_config,
    force_config,
    params=params,
    additional_dynamics=additional_dynamics,
)

# Calculate orbital period and propagate for 3 orbits
orbital_period = bh.orbital_period(oe[0])
num_orbits = 3
total_time = num_orbits * orbital_period

print(f"\nOrbital period: {orbital_period:.1f} s ({orbital_period / 60:.1f} min)")
print(f"Propagating for {num_orbits} orbits ({total_time / 60:.1f} min)")

# Propagate
prop.propagate_to(epoch + total_time)

# Check final state
final_state = prop.current_state()
final_charge = final_state[6]
charge_change = final_charge - initial_charge

print("\nFinal battery state:")
print(
    f"  Final charge: {final_charge:.2f} Wh ({100 * final_charge / battery_capacity:.1f}% SOC)"
)
print(f"  Charge change: {charge_change:+.2f} Wh")

# Sample trajectory to find eclipse statistics
traj = prop.trajectory
dt = 30.0  # 30 second samples
t = 0.0
eclipse_time = 0.0
sunlit_time = 0.0

while t <= total_time:
    current_epoch = epoch + t
    try:
        state = traj.interpolate(current_epoch)
        r_eci = state[:3]
        r_sun = bh.sun_position(current_epoch)
        illumination = bh.eclipse_conical(r_eci, r_sun)

        if illumination < 0.01:  # In eclipse (< 1% illumination)
            eclipse_time += dt
        else:
            sunlit_time += dt
    except RuntimeError:
        pass
    t += dt

eclipse_fraction = eclipse_time / (eclipse_time + sunlit_time)
print("\nEclipse statistics:")
print(
    f"  Sunlit time: {sunlit_time / 60:.1f} min ({100 * (1 - eclipse_fraction):.1f}%)"
)
print(f"  Eclipse time: {eclipse_time / 60:.1f} min ({100 * eclipse_fraction:.1f}%)")

# Validate
assert final_charge > 0, "Battery should not be depleted"
assert final_charge <= battery_capacity, "Battery should not exceed capacity"
assert eclipse_time > 0, "Should have some eclipse periods"
assert sunlit_time > 0, "Should have some sunlit periods"

print("\nExample validated successfully!")

  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
162
163
164
165
166
167
168
169
170
171
172
173
use brahe as bh;
use bh::integrators::traits::DStateDynamics;
use bh::propagators::{DNumericalOrbitPropagator, ForceModelConfig, NumericalPropagationConfig};
use bh::time::Epoch;
use bh::traits::{DStatePropagator, InterpolatableTrajectory};
use nalgebra as na;
use std::sync::Arc;

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

    // Create initial epoch
    let epoch = Epoch::from_datetime(2024, 6, 21, 12, 0, 0.0, 0.0, bh::TimeSystem::UTC);

    // Initial orbital elements and state - LEO orbit
    let oe = na::SVector::<f64, 6>::new(bh::R_EARTH + 500e3, 0.01, 45.0, 0.0, 0.0, 0.0);
    let orbital_state = bh::state_koe_to_eci(oe, bh::AngleFormat::Degrees);

    // Extended state: [x, y, z, vx, vy, vz, battery_charge]
    let battery_capacity = 100.0; // Wh
    let initial_charge = 80.0; // Wh (80% SOC)
    let mut initial_state = na::DVector::zeros(7);
    for i in 0..6 {
        initial_state[i] = orbital_state[i];
    }
    initial_state[6] = initial_charge;

    // Power system parameters
    let solar_panel_power = 50.0; // W (when fully illuminated)
    let load_power = 30.0; // W (continuous consumption)

    println!("Power system parameters:");
    println!("  Battery capacity: {} Wh", battery_capacity);
    println!(
        "  Initial charge: {} Wh ({:.0}% SOC)",
        initial_charge,
        100.0 * initial_charge / battery_capacity
    );
    println!("  Solar panel power: {} W", solar_panel_power);
    println!("  Load power: {} W", load_power);
    println!(
        "  Net charging rate (sunlit): {} W",
        solar_panel_power - load_power
    );
    println!("  Net discharge rate (eclipse): {} W", load_power);

    // Store epoch in Arc for closure
    let epoch_ref = Arc::new(epoch);
    let epoch_clone = Arc::clone(&epoch_ref);

    // Additional dynamics for battery tracking
    let additional_dynamics: DStateDynamics = Box::new(move |t, state, _params| {
        let mut dx = na::DVector::zeros(state.len());
        let r_eci = na::Vector3::new(state[0], state[1], state[2]);

        // Get sun position at current epoch
        let current_epoch = *epoch_clone + t;
        let r_sun = bh::sun_position(current_epoch);

        // Get illumination fraction (0 = umbra, 0-1 = penumbra, 1 = sunlit)
        let illumination = bh::eclipse_conical(r_eci, r_sun);

        // Battery dynamics (Wh/s = W / 3600)
        let power_in = illumination * solar_panel_power; // W
        let power_out = load_power; // W
        let mut charge_rate = (power_in - power_out) / 3600.0; // Wh/s

        // Apply battery limits (0 to capacity)
        let charge = state[6];
        if charge >= battery_capacity && charge_rate > 0.0 {
            charge_rate = 0.0; // Battery full
        } else if charge <= 0.0 && charge_rate < 0.0 {
            charge_rate = 0.0; // Battery empty
        }

        dx[6] = charge_rate;
        dx
    });

    // Create propagator with two-body dynamics
    let mut prop = DNumericalOrbitPropagator::new(
        epoch,
        initial_state.clone(),
        NumericalPropagationConfig::default(),
        ForceModelConfig::two_body_gravity(),
        None, // params
        Some(additional_dynamics),
        None, // control_input
        None, // initial_covariance
    )
    .unwrap();

    // Calculate orbital period and propagate for 3 orbits
    let orbital_period = bh::orbital_period(oe[0]);
    let num_orbits = 3;
    let total_time = num_orbits as f64 * orbital_period;

    println!(
        "\nOrbital period: {:.1} s ({:.1} min)",
        orbital_period,
        orbital_period / 60.0
    );
    println!(
        "Propagating for {} orbits ({:.1} min)",
        num_orbits,
        total_time / 60.0
    );

    // Propagate
    prop.propagate_to(epoch + total_time);

    // Check final state
    let final_state = prop.current_state();
    let final_charge = final_state[6];
    let charge_change = final_charge - initial_charge;

    println!("\nFinal battery state:");
    println!(
        "  Final charge: {:.2} Wh ({:.1}% SOC)",
        final_charge,
        100.0 * final_charge / battery_capacity
    );
    println!("  Charge change: {:+.2} Wh", charge_change);

    // Sample trajectory to find eclipse statistics
    let traj = prop.trajectory();
    let dt = 30.0; // 30 second samples
    let mut t = 0.0;
    let mut eclipse_time = 0.0;
    let mut sunlit_time = 0.0;

    while t <= total_time {
        let current_epoch = epoch + t;
        if let Ok(state) = traj.interpolate(&current_epoch) {
            let r_eci = na::Vector3::new(state[0], state[1], state[2]);
            let r_sun = bh::sun_position(current_epoch);
            let illumination = bh::eclipse_conical(r_eci, r_sun);

            if illumination < 0.01 {
                // In eclipse (< 1% illumination)
                eclipse_time += dt;
            } else {
                sunlit_time += dt;
            }
        }
        t += dt;
    }

    let eclipse_fraction = eclipse_time / (eclipse_time + sunlit_time);
    println!("\nEclipse statistics:");
    println!(
        "  Sunlit time: {:.1} min ({:.1}%)",
        sunlit_time / 60.0,
        100.0 * (1.0 - eclipse_fraction)
    );
    println!(
        "  Eclipse time: {:.1} min ({:.1}%)",
        eclipse_time / 60.0,
        100.0 * eclipse_fraction
    );

    // Validate
    assert!(final_charge > 0.0, "Battery should not be depleted");
    assert!(
        final_charge <= battery_capacity,
        "Battery should not exceed capacity"
    );
    assert!(eclipse_time > 0.0, "Should have some eclipse periods");
    assert!(sunlit_time > 0.0, "Should have some sunlit periods");

    println!("\nExample validated successfully!");
}

Battery Charge and Illumination Profile

The following plot shows battery state of charge over 3 orbits, with illumination fraction and eclipse periods clearly visible:

Plot Source
battery_tracking.py
  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
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

# Add plots directory to path for importing brahe_theme
sys.path.insert(0, str(pathlib.Path(__file__).parent.parent.parent))
from brahe_theme import save_themed_html, get_theme_colors

# Configuration
SCRIPT_NAME = pathlib.Path(__file__).stem
OUTDIR = pathlib.Path(os.getenv("BRAHE_FIGURE_OUTPUT_DIR", "./docs/figures/"))
os.makedirs(OUTDIR, exist_ok=True)

# Initialize EOP data
bh.initialize_eop()

# Create initial epoch and state - LEO orbit
epoch = bh.Epoch.from_datetime(2024, 6, 21, 12, 0, 0.0, 0.0, bh.TimeSystem.UTC)
oe = np.array([bh.R_EARTH + 500e3, 0.01, 45.0, 15.0, 30.0, 45.0])
orbital_state = bh.state_koe_to_eci(oe, bh.AngleFormat.DEGREES)

# Extended state: [x, y, z, vx, vy, vz, battery_charge]
battery_capacity = 100.0  # Wh
initial_charge = 80.0  # Wh (80% SOC)
initial_state = np.concatenate([orbital_state, [initial_charge]])

# Power system parameters
solar_panel_power = 50.0  # W (when fully illuminated)
load_power = 30.0  # W (continuous consumption)

# Spacecraft parameters for force model
params = np.array([500.0, 2.0, 2.2, 2.0, 1.3])


# Define additional dynamics for battery tracking
def additional_dynamics(t, state, params):
    dx = np.zeros(len(state))
    r_eci = state[:3]

    # Get sun position at current epoch
    current_epoch = epoch + t
    r_sun = bh.sun_position(current_epoch)

    # Get illumination fraction (0 = umbra, 0-1 = penumbra, 1 = sunlit)
    illumination = bh.eclipse_conical(r_eci, r_sun)

    # Battery dynamics (Wh/s = W / 3600)
    power_in = illumination * solar_panel_power
    power_out = load_power
    charge_rate = (power_in - power_out) / 3600.0

    # Apply battery limits
    charge = state[6]
    if charge >= battery_capacity and charge_rate > 0:
        charge_rate = 0.0
    elif charge <= 0 and charge_rate < 0:
        charge_rate = 0.0

    dx[6] = charge_rate
    return dx


# Create propagator with two-body dynamics
force_config = bh.ForceModelConfig.two_body()
prop_config = bh.NumericalPropagationConfig.default()

prop = bh.NumericalOrbitPropagator(
    epoch,
    initial_state,
    prop_config,
    force_config,
    params=params,
    additional_dynamics=additional_dynamics,
)

# Propagate for 3 orbits
orbital_period = bh.orbital_period(oe[0])
num_orbits = 3
total_time = num_orbits * orbital_period

prop.propagate_to(epoch + total_time)

# Sample trajectory for plotting
traj = prop.trajectory
times = []
charge_vals = []
illumination_vals = []

dt = 10.0  # 10 second samples
t = 0.0
while t <= total_time:
    current_epoch = epoch + t
    try:
        state = traj.interpolate(current_epoch)
        r_eci = state[:3]
        r_sun = bh.sun_position(current_epoch)
        illumination = bh.eclipse_conical(r_eci, r_sun)

        times.append(t / 60.0)  # Convert to minutes
        charge_vals.append(state[6])
        illumination_vals.append(illumination)
    except RuntimeError:
        pass

    t += dt

# Find eclipse regions for shading
eclipse_starts = []
eclipse_ends = []
in_eclipse = False
eclipse_threshold = 0.01  # Consider <1% illumination as eclipse

for i, illum in enumerate(illumination_vals):
    if illum < eclipse_threshold and not in_eclipse:
        eclipse_starts.append(times[i])
        in_eclipse = True
    elif illum >= eclipse_threshold and in_eclipse:
        eclipse_ends.append(times[i])
        in_eclipse = False

# Close last eclipse if still in eclipse at end
if in_eclipse and len(eclipse_starts) > len(eclipse_ends):
    eclipse_ends.append(times[-1])

# Calculate statistics
final_charge = charge_vals[-1]
charge_change = final_charge - initial_charge
sunlit_time = sum(1 for i in illumination_vals if i > eclipse_threshold) * dt / 60.0
eclipse_time = len(times) * dt / 60.0 - sunlit_time


def create_figure(theme):
    colors = get_theme_colors(theme)

    # Create figure with secondary y-axis
    fig = make_subplots(specs=[[{"secondary_y": True}]])

    # Battery charge (primary y-axis)
    fig.add_trace(
        go.Scatter(
            x=times,
            y=charge_vals,
            mode="lines",
            name="Battery Charge",
            line=dict(color=colors["primary"], width=3),
        ),
        secondary_y=False,
    )

    # Illumination fraction (secondary y-axis) as filled area
    fig.add_trace(
        go.Scatter(
            x=times,
            y=illumination_vals,
            mode="lines",
            name="Illumination",
            line=dict(color=colors["secondary"], width=1),
            fill="tozeroy",
            fillcolor="rgba(255, 165, 0, 0.15)"
            if theme == "light"
            else "rgba(255, 170, 68, 0.15)",
        ),
        secondary_y=True,
    )

    # Add eclipse shading
    for start, end in zip(eclipse_starts, eclipse_ends):
        fig.add_vrect(
            x0=start,
            x1=end,
            fillcolor="rgba(100, 100, 100, 0.2)"
            if theme == "light"
            else "rgba(50, 50, 50, 0.4)",
            layer="below",
            line_width=0,
        )

    # Add reference lines
    fig.add_hline(
        y=initial_charge,
        line_dash="dot",
        line_color="gray",
        annotation_text=f"Initial: {initial_charge:.0f} Wh",
        annotation_position="top right",
        secondary_y=False,
    )

    fig.add_hline(
        y=battery_capacity,
        line_dash="dot",
        line_color=colors["accent"],
        annotation_text=f"Capacity: {battery_capacity:.0f} Wh",
        annotation_position="bottom right",
        secondary_y=False,
    )

    # Update layout
    fig.update_layout(
        title=f"Battery Charge with Eclipse Cycles (LEO, {num_orbits} orbits)",
        xaxis_title="Time (min)",
        height=500,
        margin=dict(l=60, r=80, t=80, b=60),
        legend=dict(
            yanchor="top",
            y=0.99,
            xanchor="left",
            x=0.01,
            bgcolor="rgba(255,255,255,0.8)"
            if theme == "light"
            else "rgba(30,30,30,0.8)",
        ),
    )

    # Update y-axes
    fig.update_yaxes(
        title_text="Battery Charge (Wh)",
        range=[60, 105],
        secondary_y=False,
    )
    fig.update_yaxes(
        title_text="Illumination Fraction",
        range=[0, 1.1],
        secondary_y=True,
    )

    # Add summary annotation
    summary_text = (
        f"Charge change: {charge_change:+.2f} Wh<br>"
        f"Final SOC: {100 * final_charge / battery_capacity:.1f}%<br>"
        f"Eclipse: {eclipse_time:.1f} min ({100 * eclipse_time / (eclipse_time + sunlit_time):.0f}%)"
    )
    fig.add_annotation(
        x=0.98,
        y=0.02,
        xref="paper",
        yref="paper",
        text=summary_text,
        showarrow=False,
        font=dict(size=11),
        align="right",
        bordercolor="gray",
        borderwidth=1,
        borderpad=4,
        bgcolor="white" if theme == "light" else "#1e1e1e",
        opacity=0.9,
    )

    # Add "Eclipse" label to first eclipse region
    if eclipse_starts:
        mid_eclipse = (eclipse_starts[0] + eclipse_ends[0]) / 2
        fig.add_annotation(
            x=mid_eclipse,
            y=0.95,
            xref="x",
            yref="paper",
            text="Eclipse",
            showarrow=False,
            font=dict(size=10, color=colors["font_color"]),
        )

    return fig


# Save themed HTML files
light_path, dark_path = save_themed_html(create_figure, OUTDIR / SCRIPT_NAME)
print(f"Generated {light_path}")
print(f"Generated {dark_path}")

The battery charges during sunlit periods (illumination = 1) and discharges during eclipse (illumination = 0). The penumbra regions show gradual transitions in illumination.

Other Common Extensions

Attitude Dynamics

Track spacecraft attitude alongside orbital motion. This example shows quaternion attitude propagation:

 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
# State: [pos(3), vel(3), q0, q1, q2, q3, wx, wy, wz] = 13 elements
# Quaternion [q0, q1, q2, q3] + angular velocity [wx, wy, wz]
def additional_dynamics_attitude(t, state, params):
    """Attitude dynamics only - orbital handled by force model."""
    dx = np.zeros(len(state))

    # Extract quaternion and angular velocity
    q = state[6:10]  # [q0, q1, q2, q3]
    omega = state[10:13]  # [wx, wy, wz] rad/s

    # Quaternion kinematics: dq/dt = 0.5 * Omega(omega) * q
    omega_matrix = np.array([
        [0, -omega[0], -omega[1], -omega[2]],
        [omega[0], 0, omega[2], -omega[1]],
        [omega[1], -omega[2], 0, omega[0]],
        [omega[2], omega[1], -omega[0], 0]
    ])
    dq = 0.5 * omega_matrix @ q
    dx[6:10] = dq

    # Angular velocity dynamics: I * domega/dt = -omega x (I * omega) + torque
    I = np.diag([10.0, 12.0, 8.0])  # Inertia tensor (kg*m^2)
    torque = np.zeros(3)  # External torques
    domega = np.linalg.solve(I, -np.cross(omega, I @ omega) + torque)
    dx[10:13] = domega

    return dx

 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
// State: [pos(3), vel(3), q0, q1, q2, q3, wx, wy, wz] = 13 elements
let additional_dynamics_attitude: DStateDynamics = Box::new(
    |_t: f64, state: &na::DVector<f64>, _params: Option<&na::DVector<f64>>| {
        let mut dx = na::DVector::zeros(state.len());

        // Extract quaternion and angular velocity
        let q = na::Vector4::new(state[6], state[7], state[8], state[9]);
        let omega = na::Vector3::new(state[10], state[11], state[12]);

        // Quaternion kinematics: dq/dt = 0.5 * Omega(omega) * q
        let omega_matrix = na::Matrix4::new(
            0.0, -omega[0], -omega[1], -omega[2],
            omega[0], 0.0, omega[2], -omega[1],
            omega[1], -omega[2], 0.0, omega[0],
            omega[2], omega[1], -omega[0], 0.0,
        );
        let dq = 0.5 * omega_matrix * q;
        dx[6] = dq[0]; dx[7] = dq[1]; dx[8] = dq[2]; dx[9] = dq[3];

        // Angular velocity dynamics: I * domega/dt = -omega x (I * omega) + torque
        let inertia = na::Matrix3::from_diagonal(&na::Vector3::new(10.0, 12.0, 8.0));
        let torque = na::Vector3::zeros();
        let domega = inertia.try_inverse().unwrap()
            * (-omega.cross(&(inertia * omega)) + torque);
        dx[10] = domega[0]; dx[11] = domega[1]; dx[12] = domega[2];

        dx
    },
);

Thermal State

Track spacecraft temperature:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
# State: [pos(3), vel(3), temperature]
def additional_dynamics_thermal(t, state, params):
    """Thermal dynamics only - orbital handled by force model."""
    dx = np.zeros(len(state))
    temp = state[6]

    # Simplified radiation balance
    q_solar = solar_flux_absorbed(state[:3])
    q_radiated = emissivity * stefan_boltzmann * temp**4 * area
    dx[6] = (q_solar - q_radiated) / (mass * specific_heat)

    return dx

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
// State: [pos(3), vel(3), temperature]
let additional_dynamics_thermal: bh::DAdditionalDynamics = Some(Box::new(
    |_t: f64, state: &na::DVector<f64>, _params: Option<&na::DVector<f64>>| {
        let mut dx = na::DVector::zeros(state.len());
        let temp = state[6];

        // Simplified radiation balance
        let q_solar = solar_flux_absorbed(&state.fixed_rows::<3>(0).into());
        let q_radiated = emissivity * stefan_boltzmann * temp.powi(4) * area;
        dx[6] = (q_solar - q_radiated) / (mass * specific_heat);

        dx
    },
));

Multiple Extensions

State vectors can include multiple extensions:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
# State: [pos(3), vel(3), mass, battery, temperature] = 9 elements
initial_state = np.array([
    *orbital_state,  # Position and velocity (6)
    1000.0,          # Mass (kg)
    100.0,           # Battery (Wh)
    293.0,           # Temperature (K)
])

def additional_dynamics_multi(t, state, params):
    """Return full state-sized vector with derivatives for extended elements."""
    dx = np.zeros(len(state))
    mass = state[6]
    charge = state[7]
    temp = state[8]

    dx[6] = -mass_flow_rate if thrusting else 0.0
    dx[7] = solar_input - power_consumption if is_sunlit(state[:3]) else -power_consumption
    dx[8] = (q_solar - q_radiated) / (mass * specific_heat)

    return dx

 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
// State: [pos(3), vel(3), mass, battery, temperature] = 9 elements
let initial_state = na::DVector::from_vec(vec![
    orbital_state[0], orbital_state[1], orbital_state[2],  // Position
    orbital_state[3], orbital_state[4], orbital_state[5],  // Velocity
    1000.0,  // Mass (kg)
    100.0,   // Battery (Wh)
    293.0,   // Temperature (K)
]);

let additional_dynamics_multi: bh::DAdditionalDynamics = Some(Box::new(
    |_t: f64, state: &na::DVector<f64>, _params: Option<&na::DVector<f64>>| {
        let mut dx = na::DVector::zeros(state.len());
        let mass = state[6];
        let _charge = state[7];
        let _temp = state[8];

        dx[6] = if thrusting { -mass_flow_rate } else { 0.0 };
        dx[7] = if is_sunlit(&state.fixed_rows::<3>(0).into()) {
            solar_input - power_consumption
        } else {
            -power_consumption
        };
        dx[8] = (q_solar - q_radiated) / (mass * specific_heat);

        dx
    },
));

Implementation Notes

Another way to implement extended state propagation is to use NumericalPropagator, which requires implementing the full dynamics function including orbital and extended state dynamics. However, using NumericalOrbitPropagator with additional_dynamics is often more convenient for orbital applications, as it handles standard orbital perturbations automatically. See the Generic Dynamics Propagation guide for details on using NumericalPropagator which may be preferable for highly customized dynamics.

See Also