ADCS.satellite_hardware.disturbances.srp_disturbance module

class ADCS.satellite_hardware.disturbances.srp_disturbance.SRP_Disturbance(config)

Bases: Disturbance

Solar Radiation Pressure (SRP) Disturbance Torque Model

This class models the disturbance torque caused by solar radiation pressure (SRP), i.e., the momentum flux of sunlight impinging on and reflecting from spacecraft exterior surfaces.

A discrete set of planar faces is used to represent the spacecraft geometry. For each face that is illuminated (Sun-facing), a force is computed using a simplified optical reflection model with absorptive, diffuse, and specular components. The net SRP torque is the sum of moments of these per-face forces about the spacecraft center of mass (COM).

Surface geometry and optical parameters are provided by:

GeometryConfig

The Sun vector and orbital state are provided by:

Orbital_State

Reference Vectors

Let:

• \(\mathbf{S}_B\) be the Sun position/direction vector expressed in the body frame (as returned by get_state_vector())

• \(\mathbf{R}_B\) be the spacecraft position vector (Earth to spacecraft) expressed in the body frame

The implemented model forms the body-frame unit Sun-look direction:

\[\mathbf{s}_b = \frac{\mathbf{S}_B - \mathbf{R}_B}{\left\lVert \mathbf{S}_B - \mathbf{R}_B \right\rVert},\]

using normalize().

Illumination Gate

SRP is applied only if the spacecraft is illuminated. Illumination is checked via:

is_sunlit()

If the spacecraft is in eclipse (not sunlit), the SRP torque is identically zero.

Optical Coefficients

Each face \(i=1,\dots,N\) has:

• area \(A_i\) [m²]

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

• unit normal \(\mathbf{n}_i\) (body frame)

• coefficients \(\eta_{s,i}\), \(\eta_{d,i}\), \(\eta_{a,i}\)

The coefficients represent specular, diffuse, and absorptive contributions, respectively. They are typically constrained such that \(\eta_{s,i} + \eta_{d,i} + \eta_{a,i} = 1\), though this class does not enforce the constraint.

Radiation Pressure

The solar radiation pressure magnitude is modeled using the solar constant \(S_0\) and the speed of light \(c\), both provided by:

EarthConstants

\[P_\odot = \frac{S_0}{c}.\]

Implemented Force/Torque Form

Define the incidence cosine:

\[\cos\gamma_i = \mathbf{n}_i^\top \mathbf{s}_b, \qquad \cos^+\gamma_i = \max(0,\cos\gamma_i).\]

The implementation uses two scalar multipliers:

\[m_{s,i} = A_i(\eta_{a,i}+\eta_{d,i})\,\cos^+\gamma_i,\]

\[m_{n,i} = A_i\left(2\,\eta_{s,i}\,(\cos^+\gamma_i)^2 + \tfrac{2}{3}\,\eta_{d,i}\right)\cos^+\gamma_i,\]

matching the code’s use of proj_area = A_i cos^+γ_i and m_n = proj_area * (2 η_s cos^+γ_i + 2/3 η_d).

Let the lever arm from COM to the face centroid be:

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

The total SRP torque (body frame) is:

\[\mathbf{T}_{\mathrm{SRP}} = -\frac{S_0}{c}\sum_{i=1}^N \Big[ m_{s,i}\,(\mathbf{c}_i \times \mathbf{s}_b) + m_{n,i}\,(\mathbf{c}_i \times \mathbf{n}_i) \Big].\]

Differentiability

The clamping \(\cos^+\gamma_i = \max(0,\cos\gamma_i)\) introduces a kink at \(\cos\gamma_i = 0\). The Jacobian and Hessian methods provided by this class use piecewise derivatives with a Heaviside gate (subgradient-like behavior).

Parameters:

config (GeometryConfig) – Surface geometry and optical coefficients for each face.

torque(sat, x, os)

Compute the SRP torque in the body frame.

This method obtains the Sun and spacecraft position vectors in the body frame via:

get_state_vector()

using keys:

• "s" : \(\mathbf{S}_B\)

• "r" : \(\mathbf{R}_B\)

The unit Sun-look direction is:

\[\mathbf{s}_b = \frac{\mathbf{S}_B - \mathbf{R}_B}{\left\lVert \mathbf{S}_B - \mathbf{R}_B \right\rVert}.\]

Incidence for each face is computed as:

\[\cos\gamma_i = \mathbf{n}_i^\top \mathbf{s}_b, \qquad \cos^+\gamma_i = \max(0,\cos\gamma_i).\]

Using the per-face multipliers:

\[m_{s,i} = A_i(\eta_{a,i}+\eta_{d,i})\,\cos^+\gamma_i,\]

\[m_{n,i} = A_i\left(2\,\eta_{s,i}\,(\cos^+\gamma_i)^2 + \tfrac{2}{3}\,\eta_{d,i}\right)\cos^+\gamma_i,\]

and lever arms \(\mathbf{c}_i = \mathbf{r}_i - \mathbf{r}_{\mathrm{COM}}\), the torque is:

\[\mathbf{T}_{\mathrm{SRP}} = -\frac{S_0}{c}\sum_{i=1}^N \Big[ m_{s,i}\,(\mathbf{c}_i \times \mathbf{s}_b) + m_{n,i}\,(\mathbf{c}_i \times \mathbf{n}_i) \Big].\]

If the spacecraft is not illuminated (os.is_sunlit() is false), then:

\[\mathbf{T}_{\mathrm{SRP}} = \mathbf{0}.\]

Parameters:

• sat (Satellite) – Satellite instance providing the COM position via sat.COM.

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

• os (Orbital_State) – Orbital state provider supplying "s", "r", and illumination status.

Returns:

SRP disturbance torque in the body frame [N·m], shape (3,).

Return type:

numpy.ndarray

torque_qjav(sat, x, os)

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

Let \(\mathbf{q}\) denote the attitude quaternion inside \(x\). The SRP torque depends on \(\mathbf{q}\) through the body-frame Sun-look direction:

\[\mathbf{s}_b(\mathbf{q}) = \frac{\mathbf{S}_B(\mathbf{q}) - \mathbf{R}_B(\mathbf{q})} {\left\lVert \mathbf{S}_B(\mathbf{q}) - \mathbf{R}_B(\mathbf{q}) \right\rVert}.\]

The Jacobian of \(\mathbf{s}_b\) is computed using:

normed_vec_jac()

on the vector \(\mathbf{S}_B - \mathbf{R}_B\), with quaternion derivatives provided by the orbital state.

Clamped incidence derivative

For each face:

\[\cos\gamma_i(\mathbf{q}) = \mathbf{n}_i^\top \mathbf{s}_b(\mathbf{q}), \qquad \cos^+\gamma_i(\mathbf{q}) = \max(0,\cos\gamma_i(\mathbf{q})).\]

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

\[\frac{\partial \cos^+\gamma_i}{\partial \mathbf{q}} = H(\cos\gamma_i)\, \mathbf{n}_i^\top \frac{\partial \mathbf{s}_b}{\partial \mathbf{q}}.\]

Derivative of per-face multipliers

With:

\[m_{s,i} = A_i(\eta_{a,i}+\eta_{d,i})\,\cos^+\gamma_i,\]

\[m_{n,i} = A_i\left(2\,\eta_{s,i}\,(\cos^+\gamma_i)^2 + \tfrac{2}{3}\,\eta_{d,i}\right)\cos^+\gamma_i,\]

the derivatives are:

\[\frac{\partial m_{s,i}}{\partial \mathbf{q}} = A_i(\eta_{a,i}+\eta_{d,i}) \frac{\partial \cos^+\gamma_i}{\partial \mathbf{q}},\]

\[\frac{\partial m_{n,i}}{\partial \mathbf{q}} = A_i\left( 6\eta_{s,i}(\cos^+\gamma_i)^2 + \tfrac{2}{3}\eta_{d,i} \right) \frac{\partial \cos^+\gamma_i}{\partial \mathbf{q}},\]

where the factor arises from differentiating the code-equivalent scalar form \(m_{n,i} = A_i\cos^+\gamma_i\left(2\eta_{s,i}\cos^+\gamma_i + \tfrac{2}{3}\eta_{d,i}\right)\).

Torque Jacobian

With \(\mathbf{c}_i = \mathbf{r}_i - \mathbf{r}_{\mathrm{COM}}\):

\[\frac{\partial \mathbf{T}_{\mathrm{SRP}}}{\partial \mathbf{q}} = -\frac{S_0}{c}\sum_{i=1}^N \Big[ \frac{\partial m_{s,i}}{\partial \mathbf{q}}(\mathbf{c}_i\times\mathbf{s}_b) + m_{s,i}\,\mathbf{c}_i\times\frac{\partial \mathbf{s}_b}{\partial \mathbf{q}} + \frac{\partial m_{n,i}}{\partial \mathbf{q}}(\mathbf{c}_i\times\mathbf{n}_i) \Big].\]

If the spacecraft is not illuminated, the Jacobian is returned as zero.

Parameters:

• sat (Satellite) – Satellite instance providing the COM position via sat.COM.

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

• os (Orbital_State) – Orbital state provider supplying vectors and quaternion derivatives via get_state_vector() (typically "s", "r", "ds", "dr").

Returns:

Quaternion Jacobian ∂T_SRP/∂q, shape (3, 4).

Return type:

numpy.ndarray

torque_qqhess(sat, x, os)

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

This method returns the second derivative tensor:

\[\frac{\partial^2 \mathbf{T}_{\mathrm{SRP}}}{\partial \mathbf{q}^2} \in \mathbb{R}^{3\times 4\times 4}.\]

The second-order dependence arises through:

• \(\mathbf{s}_b(\mathbf{q})\) (a normalized vector)

• the clamped incidence \(\cos^+\gamma_i(\mathbf{q})\)

The second derivative of the normalized Sun-look direction is computed using:

normed_vec_hess().

Second-order structure

Let \(\mathbf{c}_i = \mathbf{r}_i - \mathbf{r}_{\mathrm{COM}}\) and define:

\[\mathbf{t}_{s,i} = \mathbf{c}_i\times\mathbf{s}_b, \qquad \mathbf{t}_{n,i} = \mathbf{c}_i\times\mathbf{n}_i.\]

Then:

\[\mathbf{T}_{\mathrm{SRP}} = -\frac{S_0}{c}\sum_{i=1}^N \bigl(m_{s,i}\mathbf{t}_{s,i} + m_{n,i}\mathbf{t}_{n,i}\bigr).\]

Differentiating twice and collecting terms yields the schematic form:

\[\frac{\partial^2 \mathbf{T}_{\mathrm{SRP}}}{\partial \mathbf{q}^2} = -\frac{S_0}{c}\sum_{i=1}^N\Big[ \frac{\partial^2 m_{s,i}}{\partial \mathbf{q}^2}\,\mathbf{t}_{s,i} + \frac{\partial m_{s,i}}{\partial \mathbf{q}}\otimes \frac{\partial \mathbf{t}_{s,i}}{\partial \mathbf{q}} + \left(\cdot\right)^\top + m_{s,i}\,\frac{\partial^2 \mathbf{t}_{s,i}}{\partial \mathbf{q}^2} + \frac{\partial^2 m_{n,i}}{\partial \mathbf{q}^2}\,\mathbf{t}_{n,i} \Big],\]

where \((\cdot)^\top\) denotes the symmetric transpose term with swapped quaternion indices.

Clamping and gating

The clamped cosine introduces a kink at \(\cos\gamma_i = 0\). The code implements a piecewise (Heaviside-gated) second derivative consistent with:

\[\frac{\partial \cos^+\gamma_i}{\partial \mathbf{q}} = H(\cos\gamma_i)\, \mathbf{n}_i^\top \frac{\partial \mathbf{s}_b}{\partial \mathbf{q}}, \qquad \frac{\partial^2 \cos^+\gamma_i}{\partial \mathbf{q}^2} \approx H(\cos\gamma_i)\, \mathbf{n}_i^\top \frac{\partial^2 \mathbf{s}_b}{\partial \mathbf{q}^2},\]

i.e., distributional terms at the kink are neglected.

If the spacecraft is not illuminated, the Hessian is returned as zero.

Parameters:

• sat (Satellite) – Satellite instance providing the COM position via sat.COM.

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

• os (Orbital_State) – Orbital state provider supplying vectors and quaternion derivatives up to second order via get_state_vector() (typically "s", "r", "ds", "dr", "dds", "ddr").

Returns:

Quaternion Hessian tensor ∂²T_SRP/∂q², shape (3, 4, 4).

Return type:

numpy.ndarray