ADCS.satellite_hardware.disturbances.drag_disturbance module

class ADCS.satellite_hardware.disturbances.drag_disturbance.Drag_Disturbance(config)

Bases: Disturbance

Aerodynamic Drag Disturbance Model

This class models the aerodynamic drag torque acting on a satellite in low Earth orbit (LEO). A discrete set of planar surface elements (faces) is used to approximate the spacecraft geometry. For each face, a drag force is computed from the relative atmospheric flow expressed in the body frame, and the net disturbance torque is obtained by summing the moments of the face forces about the spacecraft center of mass (COM).

Geometry and Parameters

The disturbance is parameterized by a configuration object:

GeometryConfig

which provides, for each face \(i = 1,\dots,N\):

• area \(A_i\) [m²]

• centroid position \(\mathbf{r}_i\) in the body frame [m]

• outward unit normal \(\mathbf{n}_i\) in the body frame (unitless)

• drag coefficient \(C_{D,i}\) (unitless)

The satellite COM is provided by:

Satellite via sat.COM.

Relative Flow

The body-frame relative velocity vector \(\mathbf{V}_b\) and density \(\rho\) are obtained from the orbital state interface:

Orbital_State

specifically via os.get_state_vector(x=...).

Face Drag Force

For a face with normal \(\mathbf{n}_i\), define the (clipped) incidence term:

\[s_i \;=\; \max\!\bigl(0,\ \mathbf{n}_i^\top \mathbf{V}_b \bigr).\]

The implemented per-face drag force is collinear with the relative flow and acts opposite the flow direction. Using the standard dynamic-pressure scaling, a compact model consistent with the implementation is:

\[\mathbf{F}_i \;=\; -\tfrac{1}{2}\,\rho\,C_{D,i}\,A_i\, s_i\,\frac{\mathbf{V}_b}{\|\mathbf{V}_b\|}.\]

(Some implementations absorb the normalization by \(\|\mathbf{V}_b\|\) into other terms. Here, \(s_i\) is linear in \(\mathbf{V}_b\) and the net torque expression below is written to match the code path.)

Net Drag Torque

Let \(\mathbf{r}_{\mathrm{COM}}\) be the spacecraft COM in body coordinates. The total drag torque about COM is:

\[\mathbf{T}_{\mathrm{drag}} \;=\; \sum_{i=1}^N \bigl(\mathbf{r}_i - \mathbf{r}_{\mathrm{COM}}\bigr)\times \mathbf{F}_i.\]

In the form used by the implementation, define

\[F_i \;=\; C_{D,i}A_i\,\max\!\bigl(0,\mathbf{n}_i^\top\mathbf{V}_b\bigr),\]

stack lever arms \(\mathbf{c}_i = \mathbf{r}_i - \mathbf{r}_{\mathrm{COM}}\), and write:

\[\mathbf{T}_{\mathrm{drag}} \;=\; -\tfrac{1}{2}\rho \left(\sum_{i=1}^N F_i\,\mathbf{c}_i\right)\times\mathbf{V}_b.\]

Differentiability

The incidence clipping \(\max(0,\cdot)\) introduces a kink at \(\mathbf{n}_i^\top\mathbf{V}_b = 0\). The Jacobian and Hessian provided by torque_qjac() and torque_qqhess() follow the piecewise derivatives using a Heaviside-gated term. At the kink, derivatives are understood in the piecewise (one-sided) sense used by the code.

Parameters:

config (GeometryConfig) – Geometry and aerodynamic coefficients for each face.

torque(sat, x, os)

Compute the aerodynamic drag torque in the body frame.

The orbital state interface is queried as:

\[\{\mathbf{V}_b,\rho\} \leftarrow \texttt{os.get\_state\_vector}(x),\]

where \(\mathbf{V}_b\) is the relative velocity in the body frame and \(\rho\) is the atmospheric density.

If the orbital-state velocity is provided in km/s, it is converted to m/s prior to use:

\[\mathbf{V}_b[\mathrm{m/s}] \;=\; 1000\,\mathbf{V}_b[\mathrm{km/s}].\]

For each face, define the clipped incidence scalar:

\[s_i \;=\; \max\!\bigl(0,\ \mathbf{n}_i^\top \mathbf{V}_b\bigr),\]

and the scalar face factor:

\[F_i \;=\; C_{D,i}A_i\,s_i.\]

Let the lever arm from COM to face centroid be:

\[\mathbf{c}_i \;=\; \mathbf{r}_i - \mathbf{r}_{\mathrm{COM}}.\]

The implemented net torque is:

\[\mathbf{T}_{\mathrm{drag}} \;=\; -\tfrac{1}{2}\rho \left(\sum_{i=1}^N F_i\,\mathbf{c}_i\right)\times\mathbf{V}_b.\]

Units

• \(\rho\) in kg/m³

• \(\mathbf{V}_b\) in m/s

• \(A_i\) in m²

• torque output in N·m

Parameters:

• sat (Satellite) – Satellite instance providing the center of mass via sat.COM.

• x (numpy.ndarray) – Full spacecraft state; must contain the attitude quaternion used by the orbital-state mapping.

• os (Orbital_State | dict\[str, numpy.ndarray\]) – Orbital state provider that supplies "v" (body-frame relative velocity) and "rho" (density) through os.get_state_vector(x=...). Some applications may pass a mapping with equivalent keys.

Returns:

Total aerodynamic drag torque in the body frame, shape (3,).

Return type:

numpy.ndarray

torque_qjac(sat, x, os)

Compute the Jacobian of drag torque with respect to the attitude quaternion.

The drag torque depends on the attitude quaternion \(\mathbf{q}\) (contained in \(x\)) through the body-frame flow \(\mathbf{V}_b(\mathbf{q})\). The orbital state provider is expected to return:

• \(\mathbf{V}_b\) via "v"

• \(\rho\) via "rho"

• \(\dfrac{\partial \mathbf{V}_b}{\partial \mathbf{q}}\) via "dv"

from:

get_state_vector().

Piecewise derivative of incidence

For each face, define:

\[s_i(\mathbf{q}) = \max\!\bigl(0,\ \mathbf{n}_i^\top \mathbf{V}_b(\mathbf{q})\bigr).\]

Using the Heaviside step function \(H(\cdot)\):

\[\frac{\partial s_i}{\partial \mathbf{q}} \;=\; H\!\bigl(\mathbf{n}_i^\top \mathbf{V}_b\bigr)\, \mathbf{n}_i^\top \frac{\partial \mathbf{V}_b}{\partial \mathbf{q}}.\]

Then:

\[F_i(\mathbf{q}) = C_{D,i}A_i\,s_i(\mathbf{q}), \qquad \frac{\partial F_i}{\partial \mathbf{q}} = C_{D,i}A_i\,\frac{\partial s_i}{\partial \mathbf{q}}.\]

Jacobian of torque

With:

\[\mathbf{C}(\mathbf{q}) = \sum_{i=1}^N F_i(\mathbf{q})\,\mathbf{c}_i,\]

the torque is:

\[\mathbf{T}_{\mathrm{drag}}(\mathbf{q}) = -\tfrac{1}{2}\rho\ \mathbf{C}(\mathbf{q}) \times \mathbf{V}_b(\mathbf{q}).\]

Applying the product rule:

\[\frac{\partial \mathbf{T}_{\mathrm{drag}}}{\partial \mathbf{q}} = -\tfrac{1}{2}\rho\left[ \left(\frac{\partial \mathbf{C}}{\partial \mathbf{q}}\right)\times\mathbf{V}_b + \mathbf{C}\times\left(\frac{\partial \mathbf{V}_b}{\partial \mathbf{q}}\right) \right].\]

The returned Jacobian has shape (3,4) corresponding to derivatives of the 3 torque components with respect to the 4 quaternion components.

Parameters:

• sat (Satellite) – Satellite instance providing the center of mass via sat.COM.

• x (numpy.ndarray) – Full spacecraft state containing the attitude quaternion.

• os (Orbital_State) – Orbital state provider supplying "v", "rho", and "dv" through get_state_vector().

Returns:

Quaternion Jacobian \(\partial \mathbf{T}_{\mathrm{drag}}/\partial \mathbf{q}\), shape (3, 4).

Return type:

numpy.ndarray

torque_qqhess(sat, x, os)

Compute the Hessian of drag torque with respect to the attitude quaternion.

This method returns the second derivative tensor of the drag torque with respect to the quaternion components, with shape (3, 4, 4). The orbital state provider is expected to return, via get_state_vector():

• \(\mathbf{V}_b\) via "v"

• \(\rho\) via "rho"

• \(\dfrac{\partial \mathbf{V}_b}{\partial \mathbf{q}}\) via "dv"

• \(\dfrac{\partial^2 \mathbf{V}_b}{\partial \mathbf{q}^2}\) via "ddv"

Face factors and their derivatives

For each face:

\[s_i(\mathbf{q}) = \max\!\bigl(0,\ \mathbf{n}_i^\top \mathbf{V}_b(\mathbf{q})\bigr), \qquad F_i(\mathbf{q}) = C_{D,i}A_i\,s_i(\mathbf{q}).\]

Using the Heaviside gate for the first derivative:

\[\frac{\partial s_i}{\partial \mathbf{q}} = H\!\bigl(\mathbf{n}_i^\top \mathbf{V}_b\bigr)\, \mathbf{n}_i^\top \frac{\partial \mathbf{V}_b}{\partial \mathbf{q}}.\]

Likewise, the second derivative is taken piecewise (gated) as implemented:

\[\frac{\partial^2 s_i}{\partial \mathbf{q}^2} \approx H\!\bigl(\mathbf{n}_i^\top \mathbf{V}_b\bigr)\, \mathbf{n}_i^\top \frac{\partial^2 \mathbf{V}_b}{\partial \mathbf{q}^2},\]

noting that distributional terms at the kink are neglected by this piecewise model.

Therefore:

\[\frac{\partial F_i}{\partial \mathbf{q}} = C_{D,i}A_i\,\frac{\partial s_i}{\partial \mathbf{q}}, \qquad \frac{\partial^2 F_i}{\partial \mathbf{q}^2} = C_{D,i}A_i\,\frac{\partial^2 s_i}{\partial \mathbf{q}^2}.\]

Hessian of torque

Define:

\[\mathbf{C}(\mathbf{q}) = \sum_{i=1}^N F_i(\mathbf{q})\,\mathbf{c}_i.\]

The torque is:

\[\mathbf{T}_{\mathrm{drag}}(\mathbf{q}) = -\tfrac{1}{2}\rho\ \mathbf{C}(\mathbf{q})\times \mathbf{V}_b(\mathbf{q}).\]

The second derivative (arranged as a tensor over quaternion components) follows from differentiating the Jacobian and grouping the mixed terms:

\[\frac{\partial^2 \mathbf{T}_{\mathrm{drag}}}{\partial \mathbf{q}^2} = -\tfrac{1}{2}\rho\Big[ \left(\frac{\partial^2 \mathbf{C}}{\partial \mathbf{q}^2}\right)\times\mathbf{V}_b + \left(\frac{\partial \mathbf{C}}{\partial \mathbf{q}}\right)\times \left(\frac{\partial \mathbf{V}_b}{\partial \mathbf{q}}\right) + \left(\left(\frac{\partial \mathbf{C}}{\partial \mathbf{q}}\right)\times \left(\frac{\partial \mathbf{V}_b}{\partial \mathbf{q}}\right)\right)^{\!\top} + \mathbf{C}\times\left(\frac{\partial^2 \mathbf{V}_b}{\partial \mathbf{q}^2}\right) \Big],\]

where \((\cdot)^{\top}\) denotes the symmetric pairing term obtained by swapping the quaternion-derivative indices (as implemented via an explicit transpose of the mixed term).

Parameters:

• sat (Satellite) – Satellite instance providing the center of mass via sat.COM.

• x (numpy.ndarray) – Full spacecraft state containing the attitude quaternion.

• os (Orbital_State) – Orbital state provider supplying "v", "rho", "dv", and "ddv" through get_state_vector().

Returns:

Quaternion Hessian tensor \(\partial^2 \mathbf{T}_{\mathrm{drag}}/\partial \mathbf{q}^2\), shape (3, 4, 4).

Return type:

numpy.ndarray