From 53a3a64122c3df794f7c4742c83a4590a11e28d3 Mon Sep 17 00:00:00 2001 From: t Date: Mon, 26 Oct 2026 15:49:59 -0400 Subject: [PATCH 1/5] Add a force interface to disturbances and Satellite.dist_forces Disturbance gains force(), force_qjac() and force_qqhess() that return zeros, mirroring the torque trio, so torque-only disturbances need no change. Satellite.dist_forces() sums the body-frame forces over the registered disturbances using the same sat/no-sat call convention as dist_torques(). The attitude dynamics are untouched; the force is for orbit propagation and inspection. Co-Authored-By: Claude Fable 5.1 --- .../disturbances/disturbance.py | 35 ++++++ .../satellite_hardware/satellite/satellite.py | 44 +++++++ .../test_disturbance_forces.py | 113 ++++++++++++++++++ 3 files changed, 192 insertions(+) create mode 100644 testing/test_disturbances/test_disturbance_forces.py diff --git a/ADCS/satellite_hardware/disturbances/disturbance.py b/ADCS/satellite_hardware/disturbances/disturbance.py index 44d84e8d..2634bfcf 100644 --- a/ADCS/satellite_hardware/disturbances/disturbance.py +++ b/ADCS/satellite_hardware/disturbances/disturbance.py @@ -135,6 +135,41 @@ def update_errors(self, j2000: float) -> None: self.last_update_time = j2000 self.update() + def force(self, *args, **kwargs) -> np.ndarray: + r"""Net force this disturbance applies to the spacecraft, in the body frame. + + Every disturbance applies a torque, and some of them (aerodynamic + drag, solar radiation pressure, thrust) also push the centre of mass. + Subclasses that model a force override this method; the base class + returns zero so a torque-only disturbance needs no change. + + The call signature mirrors :meth:`torque`: ``force(sat, x, os)`` for + disturbances that need the satellite, ``force(x, os)`` otherwise. + :meth:`~ADCS.satellite_hardware.satellite.satellite.Satellite.dist_forces` + sums these over every disturbance. + + :return: Force in body coordinates, shape ``(3,)`` [N]. + :rtype: numpy.ndarray + """ + return np.zeros(3) + + def force_qjac(self, *args, **kwargs) -> np.ndarray: + r"""Jacobian of :meth:`force` with respect to the attitude quaternion. + + Layout matches :meth:`torque_qjac`: row ``i`` is the derivative with + respect to quaternion component ``i``, shape ``(4, 3)``. Zero in the + base class. + """ + return np.zeros((4, 3)) + + def force_qqhess(self, *args, **kwargs) -> np.ndarray: + r"""Hessian of :meth:`force` with respect to the attitude quaternion. + + Layout matches :meth:`torque_qqhess`: ``(q_i, q_j, force component)``, + shape ``(4, 4, 3)``. Zero in the base class. + """ + return np.zeros((4, 4, 3)) + def parameter_process_psd(self, *, form: str = "full") -> Covariance: r"""Return the continuous PSD of estimated disturbance parameters. diff --git a/ADCS/satellite_hardware/satellite/satellite.py b/ADCS/satellite_hardware/satellite/satellite.py index fc76484f..4d1f4aba 100644 --- a/ADCS/satellite_hardware/satellite/satellite.py +++ b/ADCS/satellite_hardware/satellite/satellite.py @@ -822,6 +822,50 @@ def dist_torques(self, x: State, os: Orbital_State, dmode: Optional[ErrorMode] = return torque_total + def dist_forces(self, x: State, os: Orbital_State, dmode: Optional[ErrorMode] = None) -> np.ndarray: + r""" + Compute the total disturbance force on the spacecraft, in the body frame. + + The total is the sum over all registered disturbances: + + .. math:: + + \mathbf{F}_{dist}(\mathbf{x},\mathrm{os}) + \;=\;\sum_{i=1}^{N_d} \mathbf{F}_{dist,i}(\mathbf{x},\mathrm{os}). + + Disturbances that only produce a torque contribute zero through the + base-class :meth:`~ADCS.satellite_hardware.disturbances.disturbance.Disturbance.force`. + As with :meth:`dist_torques`, a disturbance may accept a ``sat`` keyword + argument, in which case this satellite is passed in. + + The attitude dynamics do not consume this force (it moves the centre of + mass, not the attitude); it is provided for orbit propagation and for + inspection. Divide by :attr:`mass` for an acceleration. + + :param x: State vector. + :type x: ADCS.state.State + + :param os: Orbital/environmental state. + :type os: :class:`~ADCS.orbits.orbital_state.Orbital_State` + + :param dmode: Optional error/noise mode (unused by the current force models). + :type dmode: :class:`~ADCS.satellite_hardware.errors.ErrorMode` | None + + :return: Disturbance force vector in body frame, shape ``(3,)`` (N). + :rtype: numpy.ndarray + """ + force_total = np.zeros(3) + + for j in self.disturbances: + if 'sat' in j.force.__code__.co_varnames: + f = j.force(sat=self, x=x, os=os) + else: + f = j.force(x=x, os=os) + + force_total += np.asarray(f, dtype=float).reshape(3,) + + return force_total + def act_torque(self, x: State, u: np.ndarray, os: Orbital_State, dmode: Optional[ErrorMode] = None) -> np.ndarray: r""" Compute the total actuator torque. diff --git a/testing/test_disturbances/test_disturbance_forces.py b/testing/test_disturbances/test_disturbance_forces.py new file mode 100644 index 00000000..368bb2e8 --- /dev/null +++ b/testing/test_disturbances/test_disturbance_forces.py @@ -0,0 +1,113 @@ +"""The force side of the disturbance interface. + +Every disturbance exposes ``force``, ``force_qjac`` and ``force_qqhess`` the +same way it exposes the torque trio. Torque-only disturbances inherit zeros +from the base class, and ``Satellite.dist_forces`` sums whatever the +registered disturbances return. +""" + +import numpy as np + +from ADCS.satellite_hardware.disturbances import ( + Dipole_Disturbance, + Disturbance, + GG_Disturbance, + Prop_Disturbance, + Torque_Disturbance, +) +from testing.test_disturbances._helpers import ( + make_orbital_state, + make_satellite, + make_state, +) + + +class _ConstantForce(Disturbance): + """A disturbance that pushes with a fixed body-frame force and needs ``sat``.""" + + def __init__(self, force): + super().__init__() + self._force = np.asarray(force, dtype=float) + + def torque(self, sat, x, os): + return np.zeros(3) + + def force(self, sat, x, os): + assert sat is not None + return self._force + + +class _ConstantForceNoSat(Disturbance): + """Same, but with the two-argument ``(x, os)`` signature.""" + + def __init__(self, force): + super().__init__() + self._force = np.asarray(force, dtype=float) + + def torque(self, x, os): + return np.zeros(3) + + def force(self, x, os): + return self._force + + +def test_base_class_force_and_derivatives_are_zero(): + base = Disturbance() + x = make_state() + os = make_orbital_state() + + assert np.array_equal(base.force(x=x, os=os), np.zeros(3)) + assert np.array_equal(base.force_qjac(x=x, os=os), np.zeros((4, 3))) + assert np.array_equal(base.force_qqhess(x=x, os=os), np.zeros((4, 4, 3))) + + +def test_torque_only_disturbances_inherit_zero_force(): + sat = make_satellite() + x = make_state() + os = make_orbital_state() + torque_only = [ + Prop_Disturbance(np.array([1e-6, -2e-6, 3e-6])), + Dipole_Disturbance(np.array([0.01, 0.02, -0.03])), + GG_Disturbance(), + Torque_Disturbance(np.array([1e-7, 0.0, 0.0])), + ] + for dist in torque_only: + assert np.array_equal(dist.force(sat=sat, x=x, os=os), np.zeros(3)), type(dist) + assert dist.force_qjac(sat=sat, x=x, os=os).shape == (4, 3) + assert dist.force_qqhess(sat=sat, x=x, os=os).shape == (4, 4, 3) + + +def test_dist_forces_sums_both_call_signatures(): + f1 = np.array([1.0e-5, 0.0, -2.0e-5]) + f2 = np.array([0.0, 3.0e-5, 1.0e-5]) + sat = make_satellite( + disturbances=[ + _ConstantForce(f1), + _ConstantForceNoSat(f2), + Prop_Disturbance(np.array([1e-6, 0.0, 0.0])), + ] + ) + x = make_state() + os = make_orbital_state() + + assert np.allclose(sat.dist_forces(x=x, os=os), f1 + f2) + + +def test_dist_forces_is_zero_without_force_producing_disturbances(): + sat = make_satellite(disturbances=[Prop_Disturbance(np.array([1e-6, 0.0, 0.0])), GG_Disturbance()]) + x = make_state() + os = make_orbital_state() + + assert np.array_equal(sat.dist_forces(x=x, os=os), np.zeros(3)) + + +def test_dist_forces_does_not_change_attitude_dynamics(): + sat_with = make_satellite(disturbances=[_ConstantForce(np.array([0.1, 0.2, 0.3]))]) + sat_without = make_satellite(disturbances=[]) + x = make_state(w=np.array([0.01, -0.02, 0.03])) + os = make_orbital_state() + + xdot_with = sat_with.dynamics_core(x=x, u=np.zeros(0), orbital_state=os) + xdot_without = sat_without.dynamics_core(x=x, u=np.zeros(0), orbital_state=os) + + assert np.allclose(xdot_with, xdot_without) From 08e2c0f4b86964c9cf99684f1183404d82e5006f Mon Sep 17 00:00:00 2001 From: t Date: Mon, 26 Oct 2026 15:57:21 -0400 Subject: [PATCH 2/5] Add raw (non-renormalising) quaternion finite-difference test helpers The analytic *_qjac/*_qqhess methods differentiate with respect to the four raw quaternion components. The existing fd_quat_jac/fd_quat_hess renormalise the perturbed quaternion, which drops the radial part of the derivative and only agrees when the output is small enough to hide inside the tolerance. Co-Authored-By: Claude Fable 5.1 --- testing/test_disturbances/_helpers.py | 31 +++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/testing/test_disturbances/_helpers.py b/testing/test_disturbances/_helpers.py index c655a3a9..09140063 100644 --- a/testing/test_disturbances/_helpers.py +++ b/testing/test_disturbances/_helpers.py @@ -102,3 +102,34 @@ def resolve_method(obj, *names: str): if hasattr(obj, name): return getattr(obj, name) raise AttributeError(f"{obj!r} does not define any of {names!r}") + + +def fd_quat_jac_raw(fun, x: State, eps: float = 1.0e-7) -> np.ndarray: + """Central-difference quaternion Jacobian without renormalising ``q``. + + The analytic ``*_qjac`` methods differentiate with respect to the four + raw quaternion components. ``fd_quat_jac`` renormalises each perturbed + quaternion, which removes the radial component of the derivative, so it + only agrees with the analytic result when the output is small enough to + hide inside the tolerance. + """ + jac = np.zeros((4, 3)) + for index in range(4): + delta = np.zeros(4) + delta[index] = eps + plus = fun(State(w=x.w, q=x.q + delta)) + minus = fun(State(w=x.w, q=x.q - delta)) + jac[index, :] = (plus - minus) / (2.0 * eps) + return jac + + +def fd_quat_hess_raw(fun, x: State, eps: float = 1.0e-5) -> np.ndarray: + """Central-difference quaternion Hessian without renormalising ``q``.""" + hess = np.zeros((4, 4, 3)) + for index in range(4): + delta = np.zeros(4) + delta[index] = eps + plus = fd_quat_jac_raw(fun, State(w=x.w, q=x.q + delta), eps=eps) + minus = fd_quat_jac_raw(fun, State(w=x.w, q=x.q - delta), eps=eps) + hess[index, :, :] = (plus - minus) / (2.0 * eps) + return hess From e9c5946462351dee17e6ffe54053b565abb84c7e Mon Sep 17 00:00:00 2001 From: Niclas Scheuer Date: Mon, 26 Oct 2026 15:51:45 -0400 Subject: [PATCH 3/5] Add aerodynamic drag force, Jacobian and Hessian to Drag_Disturbance Taken from 2_2/disturbance_forces_clean (#240), drag portion only. --- .../disturbances/drag_disturbance.py | 156 ++++++++++++++++++ 1 file changed, 156 insertions(+) diff --git a/ADCS/satellite_hardware/disturbances/drag_disturbance.py b/ADCS/satellite_hardware/disturbances/drag_disturbance.py index ed4a3b8e..23314f03 100644 --- a/ADCS/satellite_hardware/disturbances/drag_disturbance.py +++ b/ADCS/satellite_hardware/disturbances/drag_disturbance.py @@ -35,6 +35,13 @@ def _drag_torque_kernel(V_B, rho, COM, normals, areas, centroids, CDs): out[2] = ct * (wx*V_B[1] - wy*V_B[0]) return out + +def _drag_force_kernel(V_B, rho, normals, areas, CDs): + """Aerodynamic drag force summed over geometry faces, in body coordinates.""" + cos_alpha = np.maximum(0.0, normals @ V_B) + face_factors = CDs * areas * cos_alpha + return -0.5 * rho * np.sum(face_factors) * V_B + if TYPE_CHECKING: from ADCS.satellite_hardware.satellite.satellite import Satellite @@ -242,6 +249,155 @@ def torque(self, sat: Satellite, x: State, os: Orbital_State | Dict[str, np.ndar return _drag_torque_kernel(V_B, rho, sat.COM, self.normals, self.areas, self.centroids, self.CDs) + + def force(self, sat: Satellite, x: State, os: Orbital_State | Dict[str, np.ndarray]) -> np.ndarray: + r"""Compute the aerodynamic drag force in the body frame. + + The orbital-state interface supplies the relative atmospheric velocity + and density through ``os.get_state_vector(x=x)``: + + - ``"v"``: body-frame relative velocity in km/s + - ``"rho"``: atmospheric mass density in kg/m³ + + Velocity is converted to m/s. For each face, the clipped incidence and + force contribution are: + + .. math:: + + s_i = \max(0, \mathbf{n}_i^\top\mathbf{V}_b), + \qquad + \mathbf{F}_i = -\tfrac{1}{2}\rho C_{D,i} A_i s_i\mathbf{V}_b. + + The total force is :math:`\sum_i \mathbf{F}_i` and has units N. This + returns force rather than acceleration; divide by spacecraft mass when + an acceleration is required. + + The incidence clipping is piecewise, with the force contribution from + a face set to zero when its normal points away from the incoming flow. + + :param sat: Satellite instance (accepted for consistency with the torque API). + :type sat: :class:`~ADCS.satellite_hardware.satellite.satellite.Satellite` + :param x: Spacecraft state containing the attitude quaternion used to + express the flow in body coordinates. + :type x: :class:`~ADCS.state.State` + :param os: Orbital state providing relative velocity and atmospheric density. + :type os: :class:`~ADCS.orbits.orbital_state.Orbital_State` or dict + :return: Total aerodynamic force in body coordinates, shape ``(3,)`` [N]. + :rtype: :class:`numpy.ndarray` + """ + vecs = os.get_state_vector(x=x) + V_B = np.asarray(vecs["v"], dtype=float).reshape(3) * 1000.0 + rho = float(vecs["rho"]) + return _drag_force_kernel(V_B, rho, self.normals, self.areas, self.CDs) + + def force_qjac(self, sat: Satellite, x: State, os: Orbital_State) -> np.ndarray: + r"""Compute the Jacobian of drag force with respect to attitude quaternion. + + The orbital state supplies body-frame velocity ``"v"`` [km/s], density + ``"rho"`` [kg/m³], and its quaternion derivative ``"dv"`` through + :meth:`~ADCS.orbits.orbital_state.Orbital_State.get_state_vector`. + Velocity and its derivative are converted to m/s before evaluating drag. + + For face incidence + + .. math:: + + s_i = \max(0, \mathbf{n}_i^\top\mathbf{V}_b), + + the piecewise derivative is gated by whether the face is exposed to the + flow. Defining + + .. math:: + + C = \sum_i C_{D,i}A_i s_i, + \qquad \mathbf{F}_{drag}=-\tfrac{1}{2}\rho C\mathbf{V}_b, + + the Jacobian applies the product rule to both :math:`C` and + :math:`\mathbf{V}_b`. Incidence-boundary distribution terms are omitted. + + The returned tensor uses quaternion component first and force component + last, shape ``(4, 3)``. Its units are N per quaternion component. + + :param sat: Satellite instance (accepted for consistency with the torque API). + :type sat: :class:`~ADCS.satellite_hardware.satellite.satellite.Satellite` + :param x: Spacecraft state containing the attitude quaternion. + :type x: :class:`~ADCS.state.State` + :param os: Orbital state supplying ``"v"``, ``"rho"``, and ``"dv"``. + :type os: :class:`~ADCS.orbits.orbital_state.Orbital_State` + :return: Quaternion Jacobian :math:`\partial\mathbf{F}_{drag}/\partial\mathbf{q}`, shape ``(4, 3)``. + :rtype: :class:`numpy.ndarray` + """ + vecs = os.get_state_vector(x=x) + V_B = np.asarray(vecs["v"], dtype=float).reshape(3) * 1000.0 + rho = float(vecs["rho"]) + dV__dq = np.asarray(vecs["dv"], dtype=float).reshape(4, 3) * 1000.0 + + incidence = self.normals @ V_B + face_lit = incidence > 0.0 + weights = self.CDs * self.areas + face_factors = weights * np.maximum(incidence, 0.0) + total_factor = float(np.sum(face_factors)) + + # d(sum_i CD_i A_i max(n_i . V, 0))/dq, shape (4,). + dtotal_factor__dq = (dV__dq @ self.normals.T) @ (weights * face_lit) + return -0.5 * rho * ( + dtotal_factor__dq[:, np.newaxis] * V_B[np.newaxis, :] + + total_factor * dV__dq + ) + + def force_qqhess(self, sat: Satellite, x: State, os: Orbital_State) -> np.ndarray: + r"""Compute the Hessian of drag force with respect to attitude quaternion. + + The orbital state supplies body-frame velocity ``"v"`` [km/s], density + ``"rho"`` [kg/m³], and its first and second quaternion derivatives + ``"dv"`` and ``"ddv"``. Velocity and both derivatives are converted to + m/s, consistent with :meth:`force`. + + Let :math:`C=\sum_i C_{D,i}A_i\max(0,\mathbf{n}_i^\top\mathbf{V}_b)`. + The force is :math:`-\tfrac{1}{2}\rho C\mathbf{V}_b`; its Hessian + contains the second derivative of :math:`C`, both mixed first-derivative + terms, and the second derivative of velocity. Derivatives of the clipped + face incidence are evaluated piecewise: exposed faces contribute and + wake-side faces contribute zero. Distributional terms at zero incidence + are neglected. + + The returned tensor is ordered ``(q_i, q_j, force_component)`` with + shape ``(4, 4, 3)`` and units N per pair of quaternion components. + + :param sat: Satellite instance (accepted for consistency with the torque API). + :type sat: :class:`~ADCS.satellite_hardware.satellite.satellite.Satellite` + :param x: Spacecraft state containing the attitude quaternion. + :type x: :class:`~ADCS.state.State` + :param os: Orbital state supplying ``"v"``, ``"rho"``, ``"dv"``, and ``"ddv"``. + :type os: :class:`~ADCS.orbits.orbital_state.Orbital_State` + :return: Quaternion Hessian :math:`\partial^2\mathbf{F}_{drag}/\partial\mathbf{q}^2`, shape ``(4, 4, 3)``. + :rtype: :class:`numpy.ndarray` + """ + vecs = os.get_state_vector(x=x) + V_B = np.asarray(vecs["v"], dtype=float).reshape(3) * 1000.0 + rho = float(vecs["rho"]) + dV__dq = np.asarray(vecs["dv"], dtype=float).reshape(4, 3) * 1000.0 + ddV__dqdq = np.asarray(vecs["ddv"], dtype=float).reshape(4, 4, 3) * 1000.0 + + incidence = self.normals @ V_B + face_lit = incidence > 0.0 + weights = self.CDs * self.areas + face_factors = weights * np.maximum(incidence, 0.0) + total_factor = float(np.sum(face_factors)) + + dincidence__dq = dV__dq @ self.normals.T # (4, Nfaces) + ddincidence__dqdq = np.einsum("ija,na->ijn", ddV__dqdq, self.normals) + dtotal_factor__dq = dincidence__dq @ (weights * face_lit) + ddtotal_factor__dqdq = np.einsum( + "ijn,n->ij", ddincidence__dqdq, weights * face_lit + ) + + return -0.5 * rho * ( + ddtotal_factor__dqdq[:, :, np.newaxis] * V_B[np.newaxis, np.newaxis, :] + + dtotal_factor__dq[:, np.newaxis, np.newaxis] * dV__dq[np.newaxis, :, :] + + dtotal_factor__dq[np.newaxis, :, np.newaxis] * dV__dq[:, np.newaxis, :] + + total_factor * ddV__dqdq + ) def torque_qjac(self, sat: Satellite, x: State, os: Orbital_State) -> np.ndarray: r""" From 9669b6e01f2c2b5bbf6789cecb830bd3f776c6bb Mon Sep 17 00:00:00 2001 From: t Date: Mon, 26 Oct 2026 15:51:45 -0400 Subject: [PATCH 4/5] Test the drag force against closed forms, the torque and finite differences Adds fd_quat_jac_raw/fd_quat_hess_raw helpers that do not renormalise the perturbed quaternion; the analytic Jacobians differentiate with respect to the raw components, so the renormalising helpers only agree when the output is small enough to hide inside the tolerance. Co-Authored-By: Claude Fable 5.1 --- .../test_disturbances/drag/test_drag_force.py | 147 ++++++++++++++++++ 1 file changed, 147 insertions(+) create mode 100644 testing/test_disturbances/drag/test_drag_force.py diff --git a/testing/test_disturbances/drag/test_drag_force.py b/testing/test_disturbances/drag/test_drag_force.py new file mode 100644 index 00000000..9c528c02 --- /dev/null +++ b/testing/test_disturbances/drag/test_drag_force.py @@ -0,0 +1,147 @@ +"""Aerodynamic drag force (the torque's companion). + +The force is the per-face flat-plate drag summed without lever arms, so the +torque must equal the sum of (centroid - COM) x (per-face force). That +consistency check ties the new force to the already-tested torque. +""" + +import numpy as np + +from ADCS.satellite_hardware.disturbances import Drag_Disturbance, GeometryFace +from testing.test_disturbances._helpers import ( + fd_quat_hess_raw, + fd_quat_jac_raw, + make_geometry_config, + make_orbital_state, + make_satellite, + make_state, +) + + +def per_face_drag_forces(drag: Drag_Disturbance, x, os) -> np.ndarray: + vecs = os.get_state_vector(x=x) + velocity = vecs["v"] * 1000.0 + rho = vecs["rho"] + face_factors = np.maximum(0.0, drag.normals @ velocity) * drag.areas * drag.CDs + return -0.5 * rho * face_factors[:, None] * velocity[None, :] + + +def test_drag_force_single_face_matches_closed_form(): + drag = Drag_Disturbance( + make_geometry_config( + GeometryFace(area=2.5, centroid=np.array([0.1, 0.3, -0.2]), normal=np.array([0.0, 1.0, 0.0]), CD=1.7) + ) + ) + sat = make_satellite(disturbances=[drag]) + x = make_state() + os = make_orbital_state(V=np.array([0.0, 7.4, 0.0]), rho=4.0e-12) + + force = drag.force(sat=sat, x=x, os=os) + + speed = 7400.0 + expected = -0.5 * 4.0e-12 * 1.7 * 2.5 * speed**2 * np.array([0.0, 1.0, 0.0]) + assert np.allclose(force, expected) + assert force @ np.array([0.0, 1.0, 0.0]) < 0.0 # opposes the velocity + + +def test_drag_force_zero_density_is_zero(): + drag = Drag_Disturbance( + make_geometry_config( + GeometryFace(area=2.0, centroid=np.zeros(3), normal=np.array([0.0, 1.0, 0.0]), CD=1.4) + ) + ) + sat = make_satellite(disturbances=[drag]) + os = make_orbital_state(V=np.array([0.2, 7.5, -0.1]), rho=0.0) + + assert np.allclose(drag.force(sat=sat, x=make_state(), os=os), np.zeros(3)) + + +def test_drag_force_back_facing_face_contributes_zero(): + drag = Drag_Disturbance( + make_geometry_config( + GeometryFace(area=2.0, centroid=np.zeros(3), normal=np.array([0.0, -1.0, 0.0]), CD=1.4) + ) + ) + sat = make_satellite(disturbances=[drag]) + os = make_orbital_state(V=np.array([0.1, 7.5, 0.2]), rho=3.0e-12) + + assert np.allclose(drag.force(sat=sat, x=make_state(), os=os), np.zeros(3)) + + +def test_drag_force_opposing_faces_only_count_the_ram_face(): + ram = GeometryFace(area=2.0, centroid=np.array([0.0, 0.2, 0.0]), normal=np.array([0.0, 1.0, 0.0]), CD=2.2) + wake = GeometryFace(area=2.0, centroid=np.array([0.0, -0.2, 0.0]), normal=np.array([0.0, -1.0, 0.0]), CD=2.2) + both = Drag_Disturbance(make_geometry_config(ram, wake)) + ram_only = Drag_Disturbance(make_geometry_config(ram)) + sat = make_satellite(disturbances=[both]) + x = make_state(q=np.array([0.9, 0.1, -0.2, 0.3])) + os = make_orbital_state(V=np.array([0.3, 7.4, -0.5]), rho=4.0e-12) + + assert np.allclose(both.force(sat=sat, x=x, os=os), ram_only.force(sat=sat, x=x, os=os)) + + +def test_drag_torque_is_lever_arm_cross_per_face_force(): + drag = Drag_Disturbance( + make_geometry_config( + GeometryFace(area=2.5, centroid=np.array([0.1, 0.3, -0.2]), normal=np.array([0.0, 1.0, 0.0]), CD=1.7), + GeometryFace(area=1.0, centroid=np.array([-0.3, 0.0, 0.1]), normal=np.array([1.0, 1.0, 0.0]), CD=2.0), + ) + ) + sat = make_satellite(COM=np.array([0.02, -0.01, 0.04]), disturbances=[drag]) + x = make_state(q=np.array([0.7, -0.2, 0.3, 0.1])) + os = make_orbital_state(V=np.array([0.2, 7.4, -0.1]), rho=4.0e-12) + + per_face = per_face_drag_forces(drag, x, os) + expected_torque = np.sum(np.cross(drag.centroids - sat.COM, per_face), axis=0) + + assert np.allclose(drag.force(sat=sat, x=x, os=os), per_face.sum(axis=0)) + assert np.allclose(drag.torque(sat=sat, x=x, os=os), expected_torque) + + +def test_drag_force_qjac_matches_finite_difference(): + drag = Drag_Disturbance( + make_geometry_config( + GeometryFace(area=2.0, centroid=np.array([0.0, 0.3, 0.0]), normal=np.array([0.0, 2.0, 0.0]), CD=1.4), + GeometryFace(area=0.5, centroid=np.array([0.2, 0.0, 0.0]), normal=np.array([1.0, 0.0, 0.5]), CD=2.0), + ) + ) + sat = make_satellite(disturbances=[drag]) + x = make_state(q=np.array([0.6, 0.2, -0.1, 0.4])) + os = make_orbital_state(V=np.array([0.01, 7.5, 0.2]), rho=3.0e-12) + expected = fd_quat_jac_raw(lambda xx: drag.force(sat=sat, x=xx, os=os), x) + + jacobian = drag.force_qjac(sat=sat, x=x, os=os) + + assert jacobian.shape == (4, 3) + assert np.allclose(jacobian, expected, atol=1.0e-9, rtol=1.0e-6) + + +def test_drag_force_qqhess_matches_finite_difference(): + drag = Drag_Disturbance( + make_geometry_config( + GeometryFace(area=2.0, centroid=np.array([0.0, 0.3, 0.0]), normal=np.array([0.0, 2.0, 0.0]), CD=1.4), + GeometryFace(area=0.5, centroid=np.array([0.2, 0.0, 0.0]), normal=np.array([1.0, 0.0, 0.5]), CD=2.0), + ) + ) + sat = make_satellite(disturbances=[drag]) + x = make_state(q=np.array([0.6, 0.2, -0.1, 0.4])) + os = make_orbital_state(V=np.array([0.01, 7.5, 0.2]), rho=3.0e-12) + expected = fd_quat_hess_raw(lambda xx: drag.force(sat=sat, x=xx, os=os), x) + + hessian = drag.force_qqhess(sat=sat, x=x, os=os) + + assert hessian.shape == (4, 4, 3) + assert np.allclose(hessian, expected, atol=1.0e-7, rtol=1.0e-5) + + +def test_drag_force_reaches_satellite_dist_forces(): + drag = Drag_Disturbance( + make_geometry_config( + GeometryFace(area=2.5, centroid=np.array([0.1, 0.3, -0.2]), normal=np.array([0.0, 1.0, 0.0]), CD=1.7) + ) + ) + sat = make_satellite(disturbances=[drag]) + x = make_state() + os = make_orbital_state(V=np.array([0.2, 7.4, -0.1]), rho=4.0e-12) + + assert np.allclose(sat.dist_forces(x=x, os=os), drag.force(sat=sat, x=x, os=os)) From 726a7a08a21a5fde005717424ee9bc1a6c7f650b Mon Sep 17 00:00:00 2001 From: nscheuer Date: Mon, 26 Oct 2026 22:53:45 +0100 Subject: [PATCH 5/5] fixing [P2] The force method advertises dictionary input but cannot handle it. [drag_disturbance.py (line 253)](/home/nicla/Generalized_ADCS/ADCS/satellite_hardware/disturbances/drag_disturbance.py:253) annotates os as Orbital_State | Dict[...], but [line 288 (line 288)](/home/nicla/Generalized_ADCS/ADCS/satellite_hardware/disturbances/drag_disturbance.py:288) unconditionally calls os.get_state_vector(...). Passing a dictionary as the annotation permits raises AttributeError. Either support a mapping here or remove it from the annotation and docstring. The existing torque method has the same mismatch. [P3] The raw finite-difference helper drops other state fields. [\\_helpers.py (line 120)](/home/nicla/Generalized_ADCS/testing/test_disturbances/_helpers.py:120) rebuilds a plain State with only w and q. That discards wheel momentum and, for EstimatorState, bias and disturbance-parameter values. A function that depends on those fields will be evaluated at a different state during the finite difference. Copy x and perturb its quaternion instead. --- .../disturbances/drag_disturbance.py | 4 ++-- testing/test_disturbances/_helpers.py | 16 ++++++++++++---- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/ADCS/satellite_hardware/disturbances/drag_disturbance.py b/ADCS/satellite_hardware/disturbances/drag_disturbance.py index 23314f03..840be0b3 100644 --- a/ADCS/satellite_hardware/disturbances/drag_disturbance.py +++ b/ADCS/satellite_hardware/disturbances/drag_disturbance.py @@ -250,7 +250,7 @@ def torque(self, sat: Satellite, x: State, os: Orbital_State | Dict[str, np.ndar return _drag_torque_kernel(V_B, rho, sat.COM, self.normals, self.areas, self.centroids, self.CDs) - def force(self, sat: Satellite, x: State, os: Orbital_State | Dict[str, np.ndarray]) -> np.ndarray: + def force(self, sat: Satellite, x: State, os: Orbital_State) -> np.ndarray: r"""Compute the aerodynamic drag force in the body frame. The orbital-state interface supplies the relative atmospheric velocity @@ -281,7 +281,7 @@ def force(self, sat: Satellite, x: State, os: Orbital_State | Dict[str, np.ndarr express the flow in body coordinates. :type x: :class:`~ADCS.state.State` :param os: Orbital state providing relative velocity and atmospheric density. - :type os: :class:`~ADCS.orbits.orbital_state.Orbital_State` or dict + :type os: :class:`~ADCS.orbits.orbital_state.Orbital_State` :return: Total aerodynamic force in body coordinates, shape ``(3,)`` [N]. :rtype: :class:`numpy.ndarray` """ diff --git a/testing/test_disturbances/_helpers.py b/testing/test_disturbances/_helpers.py index 09140063..d21d3899 100644 --- a/testing/test_disturbances/_helpers.py +++ b/testing/test_disturbances/_helpers.py @@ -117,8 +117,12 @@ def fd_quat_jac_raw(fun, x: State, eps: float = 1.0e-7) -> np.ndarray: for index in range(4): delta = np.zeros(4) delta[index] = eps - plus = fun(State(w=x.w, q=x.q + delta)) - minus = fun(State(w=x.w, q=x.q - delta)) + x_plus = x.copy() + x_minus = x.copy() + x_plus.q = x.q + delta + x_minus.q = x.q - delta + plus = fun(x_plus) + minus = fun(x_minus) jac[index, :] = (plus - minus) / (2.0 * eps) return jac @@ -129,7 +133,11 @@ def fd_quat_hess_raw(fun, x: State, eps: float = 1.0e-5) -> np.ndarray: for index in range(4): delta = np.zeros(4) delta[index] = eps - plus = fd_quat_jac_raw(fun, State(w=x.w, q=x.q + delta), eps=eps) - minus = fd_quat_jac_raw(fun, State(w=x.w, q=x.q - delta), eps=eps) + x_plus = x.copy() + x_minus = x.copy() + x_plus.q = x.q + delta + x_minus.q = x.q - delta + plus = fd_quat_jac_raw(fun, x_plus, eps=eps) + minus = fd_quat_jac_raw(fun, x_minus, eps=eps) hess[index, :, :] = (plus - minus) / (2.0 * eps) return hess