satellite

class Satellite

#include <satellite.h>

Central spacecraft model integrating dynamics, actuators, and disturbances.

The Satellite class encapsulates the complete spacecraft model including:

• Attitude dynamics (quaternion and angular velocity)

• Reaction wheel momentum states

• Actuator (magnetorquer and reaction wheel) models

• Disturbance models (drag, SRP, gravity gradient)

• Cost and constraint functions for trajectory optimization

State vector structure:

\[\begin{split} \mathbf{x} = \begin{bmatrix} \boldsymbol{\omega} \\ \mathbf{q} \\ \mathbf{h} \end{bmatrix} \quad \text{where} \quad \boldsymbol{\omega} = \text{angular velocity (3D)}, \mathbf{q} = \text{attitude quaternion (4D)}, \mathbf{h} = \text{reaction wheel momentum} \end{split}\]

Control input structure:

\[\begin{split} \mathbf{u} = \begin{bmatrix} \mathbf{m} \\ \boldsymbol{\tau}_{rw} \end{bmatrix} \end{split}\]

where \(\mathbf{m}\) are magnetorquer dipoles and \(\boldsymbol{\tau}_{rw}\) are RW torques.

Public Types

using Vec3 = Eigen::Vector3d

using Vec4 = Eigen::Vector4d

using Vec7 = Eigen::Matrix<double, 7, 1>

using Mat13 = Eigen::Matrix<double, 1, 3>

using Mat33 = Eigen::Matrix3d

using Mat34 = Eigen::Matrix<double, 3, 4>

using Mat43 = Eigen::Matrix<double, 4, 3>

using Mat73 = Eigen::Matrix<double, 7, 3>

using Mat77 = Eigen::Matrix<double, 7, 7>

using VecX = Eigen::Matrix<double, Eigen::Dynamic, 1, 0, MaxDim, 1>

using MatX = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, MaxDim, MaxDim>

using DynHessXX = saltro::math::Tensor3<saltro::limits::MAX_STATE_DIM, saltro::limits::MAX_STATE_DIM, saltro::limits::MAX_STATE_DIM>

using DynHessUX = saltro::math::Tensor3<saltro::limits::MAX_CTRL_DIM, saltro::limits::MAX_STATE_DIM, saltro::limits::MAX_STATE_DIM>

using DynHessUU = saltro::math::Tensor3<saltro::limits::MAX_CTRL_DIM, saltro::limits::MAX_CTRL_DIM, saltro::limits::MAX_STATE_DIM>

using ConstrHessXX = saltro::math::Tensor3<saltro::limits::MAX_STATE_DIM, saltro::limits::MAX_STATE_DIM, saltro::limits::MAX_CONSTRAINT_DIM>

using ConstrHessUX = saltro::math::Tensor3<saltro::limits::MAX_CTRL_DIM, saltro::limits::MAX_STATE_DIM, saltro::limits::MAX_CONSTRAINT_DIM>

using ConstrHessUU = saltro::math::Tensor3<saltro::limits::MAX_CTRL_DIM, saltro::limits::MAX_CTRL_DIM, saltro::limits::MAX_CONSTRAINT_DIM>

Public Functions

Satellite()

Default constructor; creates zero inertia and empty actuators.

Satellite(const Mat33 &Jcom_in, const PlannerSettings &settings)

Construct a satellite with specified inertia and planner settings.

Parameters:

• Jcom_in – 3×3 inertia matrix (center of mass, body frame).

• settings – Planner settings (cost/constraint configurations).

~Satellite() = default

void setInertia(const Mat33 &Jcom_in)

Set the satellite’s moment of inertia.

Updates inertia and recomputes inertia without reaction wheels if needed.

Parameters:

Jcom_in – 3×3 inertia matrix.

inline const Mat33 &inertia() const

Get total inertia (including reaction wheel contribution).

Returns:

Const reference to inertia matrix.

inline const Mat33 &invInertia() const

Get inverse of total inertia.

Returns:

Const reference to inverse inertia matrix.

inline const Mat33 &inertiaNoRW() const

Get inertia contribution excluding reaction wheels.

Returns:

Const reference to spacecraft-body inertia.

inline const Mat33 &invInertiaNoRW() const

Get inverse inertia excluding reaction wheels.

Returns:

Const reference to inverse inertia (no RW).

void setGeometryConfig(const saltro::disturbances::GeometryConfig &config)

Set satellite geometry configuration.

Used by disturbance models (drag, SRP) for surface calculations.

Parameters:

config – Geometry configuration with faces.

inline const saltro::disturbances::GeometryConfig &geometryConfig() const

Get geometry configuration (const).

Returns:

Const reference to geometry.

inline saltro::disturbances::GeometryConfig &geometryConfig()

Get geometry configuration (mutable).

Returns:

Reference to geometry.

void addMTQ(const Vec3 &axis, double max_dipole)

Add a magnetorquer to the satellite.

Parameters:

• axis – Unit vector specifying coil direction (body frame).

• max_dipole – Maximum dipole moment (A·m²).

Throws:

std::runtime_error – If max MTQ count exceeded.

void addRW(const Vec3 &axis, double max_torque, double J, double h0, double h_max)

Add a reaction wheel to the satellite.

Parameters:

• axis – Unit vector specifying wheel spin axis (body frame).

• max_torque – Maximum torque output (N·m).

• J – Wheel moment of inertia (kg·m²).

• h0 – Initial momentum (N·m·s).

• h_max – Maximum momentum (N·m·s).

Throws:

std::runtime_error – If max RW count exceeded.

inline int numMTQ() const

Get number of magnetorquers.

Returns:

Number of MTQs.

inline int numRW() const

Get number of reaction wheels.

Returns:

Number of RWs.

inline int controlDim() const

Get total control input dimension.

Equals \(\text{numMTQ} + \text{numRW}\).

Returns:

Control dimension.

inline int stateDim() const

Get full state dimension.

Equals \(7 + \text{numRW}\) (AV + Q + RW momenta).

Returns:

State dimension.

inline int reducedStateDim() const

Get reduced state dimension for optimization.

Equals \(6 + \text{numRW}\) (without full attitude).

Returns:

Reduced state dimension.

const MTQ &getMTQ(int i) const

Get a magnetorquer by index (const).

Parameters:

i – MTQ index.

Returns:

Const reference to MTQ.

const RW &getRW(int i) const

Get a reaction wheel by index (const).

Parameters:

i – RW index.

Returns:

Const reference to RW.

MTQ &getMTQ(int i)

Get a magnetorquer by index (mutable).

Parameters:

i – MTQ index.

Returns:

Reference to MTQ.

RW &getRW(int i)

Get a reaction wheel by index (mutable).

Parameters:

i – RW index.

Returns:

Reference to RW.

inline void setSettings(const PlannerSettings &settings)

Set planner settings.

Parameters:

settings – Settings object.

inline const PlannerSettings &settings() const

Get planner settings (const).

Returns:

Const reference to settings.

Vec3 actuatorTorque(const VecX &x, const VecX &u, const Vec3 &B_eci) const

Compute total actuator torque.

Sums torques from all MTQs and RWs given control input and state.

Parameters:

• x – State vector (contains quaternion for ECI→body transformation).

• u – Control input.

• B_eci – Magnetic field in ECI frame (Tesla).

Returns:

3D actuator torque in body frame.

Vec3 disturbanceTorque(const VecX &x, const DisturbanceConfig &dist, const Vec3 &R_eci, const Vec3 &B_eci, const Vec3 &S_eci, const Vec3 &V_eci, const int rho) const

Compute total disturbance torque.

Sums torques from all enabled disturbance models.

Parameters:

• x – State vector.

• dist – Disturbance configuration.

• R_eci – Position in ECI frame (meters).

• B_eci – Magnetic field in ECI (Tesla).

• S_eci – Spacecraft-to-Sun in ECI (meters).

• V_eci – Velocity in ECI (m/s).

• rho – Atmospheric density (kg/m³).

Returns:

3D disturbance torque in body frame.

VecX dynamics(const VecX &x, const VecX &u, const DisturbanceConfig &dist, const Vec3 &R_eci, const Vec3 &B_eci, const Vec3 &S_eci, const Vec3 &V_eci, const int rho) const

Compute spacecraft dynamics.

Returns state derivative:

\[ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, \mathbf{u}, \text{env}) \]

Parameters:

• x – State vector.

• u – Control input.

• dist – Disturbance configuration.

• R_eci – Orbital position (ECI).

• B_eci – Magnetic field (ECI).

• S_eci – Sun direction (ECI).

• V_eci – Velocity (ECI).

• rho – Atmospheric density.

Returns:

State derivative vector.

std::tuple<MatX, MatX, MatX> dynamicsJacobians(const VecX &x, const VecX &u, const DisturbanceConfig &dist, const Vec3 &R_eci, const Vec3 &B_eci, const Vec3 &S_eci, const Vec3 &V_eci) const

Compute dynamics Jacobians.

Returns \(\frac{\partial \mathbf{f}}{\partial \mathbf{x}}\), \(\frac{\partial \mathbf{f}}{\partial \mathbf{u}}\), and \(\frac{\partial \mathbf{f}}{\partial \text{disturbances}}\).

Parameters:

• x – State.

• u – Control.

• dist – Disturbances.

• R_eci – Position vector (ECI).

• B_eci – Magnetic field (ECI).

• S_eci – Sun direction (ECI).

• V_eci – Velocity (ECI).

Returns:

Tuple of Jacobian matrices (Fx, Fu, Fd).

std::tuple<DynHessXX, DynHessUX, DynHessUU> dynamicsHessians(const VecX &x, const VecX &u, const DisturbanceConfig &dist, const Vec3 &R_eci, const Vec3 &B_eci, const Vec3 &S_eci, const Vec3 &V_eci) const

Compute dynamics Hessians (second derivatives).

Parameters:

• x – State.

• u – Control.

• dist – Disturbances.

• R_eci – Position vector (ECI).

• B_eci – Magnetic field (ECI).

• S_eci – Sun direction (ECI).

• V_eci – Velocity (ECI).

Returns:

Tuple of Hessian tensors (Hxx, Hux, Huu).

double stageCost(int k, int N, const VecX &x, const VecX &u, const Vec3 &boresight_body, const Vec4 &attitude_target, const Vec3 &B_eci, const CostConfig &cost_cfg) const

Compute stage cost (intermediate time step).

Evaluates the cost function at time step k:

\[ J_k = c(\mathbf{x}_k, \mathbf{u}_k, k) \]

Parameters:

• k – Current time step index.

• N – Total time steps.

• x – State at step k.

• u – Control at step k.

• boresight_body – Boresight direction in body frame.

• attitude_target – Attitude target (quaternion or [NaN, ECI direction]).

• B_eci – Magnetic field (ECI).

• cost_cfg – Cost configuration.

Returns:

Scalar cost value.

double terminalCost(const VecX &x, const Vec3 &boresight_body, const Vec4 &attitude_target, const Vec3 &B_eci, const CostConfig &cost_cfg) const

Compute terminal cost (final time step).

Parameters:

• x – Final state.

• boresight_body – Boresight direction in body frame.

• attitude_target – Attitude target (quaternion or [NaN, ECI direction]).

• B_eci – Magnetic field.

• cost_cfg – Cost configuration.

Returns:

Scalar terminal cost.

std::tuple<VecX, MatX, MatX> stageCostJacobians(int k, int N, const VecX &x, const VecX &u, const Vec3 &boresight_body, const Vec4 &attitude_target, const Vec3 &B_eci, const CostConfig &cost_cfg) const

Compute stage cost Jacobians.

Returns gradients of cost:

\[\frac{\partial J_k}{\partial \mathbf{x}}, \quad \frac{\partial J_k}{\partial \mathbf{u}}\]

Parameters:

• k – Time step.

• N – Total steps.

• x – State.

• u – Control.

• boresight_body – Boresight direction in body frame.

• attitude_target – Attitude target (quaternion or [NaN, ECI direction]).

• B_eci – Magnetic field.

• cost_cfg – Cost configuration.

Returns:

Tuple of (∇_x J, ∇_u J).

std::tuple<VecX, MatX, MatX> terminalCostJacobians(const VecX &x, const Vec3 &boresight_body, const Vec4 &attitude_target, const Vec3 &B_eci, const CostConfig &cost_cfg) const

Compute terminal cost Jacobians.

Delegates to stageCostJacobians with k=0, N=1 to apply terminal weights.

Parameters:

• x – State.

• boresight_body – Boresight direction in body frame.

• attitude_target – Attitude target (quaternion or [NaN, ECI direction]).

• B_eci – Magnetic field.

• cost_cfg – Cost configuration.

Returns:

Tuple of (∇_x J, ∇_u J).

std::tuple<MatX, MatX, MatX> stageCostHessians(int k, int N, const VecX &x, const VecX &u, const Vec3 &boresight_body, const Vec4 &attitude_target, const Vec3 &B_eci, const CostConfig &cost_cfg) const

Compute stage cost Hessians.

Parameters:

• k – Time step.

• N – Total steps.

• x – State.

• u – Control.

• boresight_body – Boresight direction in body frame.

• attitude_target – Attitude target (quaternion or [NaN, ECI direction]).

• B_eci – Magnetic field.

• cost_cfg – Cost configuration.

Returns:

Tuple of Hessian matrices (Hxx, Huu, Hxu).

std::tuple<MatX, MatX, MatX> terminalCostHessians(const VecX &x, const Vec3 &boresight_body, const Vec4 &attitude_target, const Vec3 &B_eci, const CostConfig &cost_cfg) const

Compute terminal cost Hessians.

Delegates to stageCostHessians with k=0, N=1 to apply terminal weights.

Parameters:

• x – State.

• boresight_body – Boresight direction in body frame.

• attitude_target – Attitude target (quaternion or [NaN, ECI direction]).

• B_eci – Magnetic field.

• cost_cfg – Cost configuration.

Returns:

Tuple of Hessian matrices (Hxx, Huu, Hxu).

double totalCost(const Eigen::Ref<const Eigen::MatrixXd> &X, const Eigen::Ref<const Eigen::MatrixXd> &U, const Eigen::Ref<const Eigen::MatrixXd> &B, const Eigen::Ref<const Eigen::MatrixXd> &boresight, const Eigen::Ref<const Eigen::MatrixXd> &attitude_target, const CostConfig &cost_cfg) const

Compute total trajectory cost.

Sums stage costs for all intermediate steps and terminal cost at final step.

Parameters:

• X – State trajectory matrix (N × state_dim). Each row is state at step k.

• U – Control trajectory matrix ((N-1) × control_dim). Each row k is control from k to k+1.

• B – Magnetic field at each step (3 × N).

• boresight – Boresight history in body frame (3 × N), column k used at step k.

• attitude_target – Attitude target trajectory (4 × N), column k used at step k.

• cost_cfg – Cost configuration.

Returns:

Total trajectory cost J = Σ_k c(x_k, u_k) + c_terminal(x_N).

VecX constraints(int k, int N, const VecX &x, const VecX &u, const Vec3 &sun_eci, const ConstraintConfig &cnst_cfg) const

Evaluate constraints.

Returns constraint violation vector. Zero means satisfied. Constraints include angular velocity limits, actuator saturations, sun avoidance.

Constraint ordering: 1) Angular velocity magnitude constraint 2) Sun-avoidance constraint (body +X boresight) 3) If k < N-1: control and RW constraints in this order:

• MTQ upper/lower bounds (2 per MTQ)

• RW torque upper/lower bounds (2 per RW)

• RW momentum upper/lower bounds (2 per RW)

• RW stiction proxy (1 per RW)

All constraints are formulated as c <= 0.

Parameters:

• k – Time step.

• N – Total steps.

• x – State.

• u – Control.

• sun_eci – Sun direction in ECI.

• cnst_cfg – Constraint configuration.

Returns:

Constraint violation vector.

std::tuple<MatX, MatX> constraintJacobians(int k, int N, const VecX &x, const VecX &u, const Vec3 &sun_eci, const ConstraintConfig &cnst_cfg) const

Compute constraint Jacobians.

Parameters:

• k – Time step.

• N – Total steps.

• x – State.

• u – Control.

• sun_eci – Sun in ECI.

• cnst_cfg – Constraint configuration.

Returns:

Tuple of (∂c/∂x, ∂c/∂u).

std::tuple<ConstrHessUU, ConstrHessUX, ConstrHessXX> constraintHessians(int k, int N, const VecX &x, const VecX &u, const Vec3 &sun_eci, const ConstraintConfig &cnst_cfg) const

Compute constraint Hessians.

Parameters:

• k – Time step.

• N – Total steps.

• x – State.

• u – Control.

• sun_eci – Sun in ECI.

• cnst_cfg – Constraint configuration.

Returns:

Tuple of Hessian tensors (Huu, Hux, Hxx).

Public Static Attributes

static constexpr int MaxDim = saltro::limits::MAX_CONSTRAINT_DIM

static constexpr int AV_INDEX = 0

Index of angular velocity in state vector.

Angular velocity \(\boldsymbol{\omega}\) occupies indices 0–2.

static constexpr int QUAT_INDEX = 3

Index of quaternion in state vector.

Quaternion \(\mathbf{q}\) occupies indices 3–6.

static constexpr int RW_MOMENTUM_INDEX = 7

Index of first reaction wheel momentum in state vector.

Reaction wheel momenta start at index 7 (for each added RW).

Private Functions

void updateInertiaNoRW()

Recompute inertia with and without RW contributions.

std::pair<Vec4, bool> processAttitudeTarget(const Vec4 &attitude_target, const Vec3 &boresight_body, const Vec4 &q_current) const

Helper to handle dual format of attitude_target parameter.

attitude_target can be in two formats:

1. ECI goal vector: [nan, x, y, z] - position vector with NaN in first element

2. Quaternion goal vector: [q0, qx, qy, qz] - valid quaternion

Parameters:

• attitude_target – The target vector (4D)

• boresight_body – Boresight direction in body frame (used only for ECI vector format)

• q_current – Current quaternion state

Returns:

Pair<q_goal_quaternion, is_eci_format>

Private Members

Mat33 Jcom_

Mat33 invJcom_

Mat33 Jcom_noRW_

Mat33 invJcom_noRW_

std::array<std::unique_ptr<MTQ>, saltro::limits::MAX_NUM_MTQ> mtq_actuators_

std::array<std::unique_ptr<RW>, saltro::limits::MAX_NUM_RW> rw_actuators_

int num_mtq_ = 0

int num_rw_ = 0

PlannerSettings settings_

saltro::disturbances::GeometryConfig geometry_config_