Applying first suggestions.

paseos/actors/actor_builder.py

@@ -14,7 +14,7 @@
 from ..power.power_device_type import PowerDeviceType
 from ..radiation.radiation_model import RadiationModel
 from .spacecraft_body_model import SpacecraftBodyModel
-from ..attitude.attitude_model import AttitudeModel
+from ..attitude.attitude_model import AttitudeModel, TorqueDisturbanceModel
 
 
 class ActorBuilder:
@@ -565,11 +565,11 @@ def set_attitude_disturbances(
         disturbance_list = []
 
         if aerodynamic:
-            disturbance_list.append("aerodynamic")
+            disturbance_list.append(TorqueDisturbanceModel.Aerodynamic)
         if gravitational:
-            disturbance_list.append("gravitational")
+            disturbance_list.append(TorqueDisturbanceModel.Gravitational)
         if magnetic:
-            disturbance_list.append("magnetic")
+            disturbance_list.append(TorqueDisturbanceModel.Magnetic)
         # Set attitude models.
         actor._attitude_model._disturbances = disturbance_list
 
@@ -579,7 +579,7 @@ def set_attitude_model(
         actor_initial_attitude_in_rad: list[float] = [0.0, 0.0, 0.0],
         actor_initial_angular_velocity: list[float] = [0.0, 0.0, 0.0],
         actor_pointing_vector_body: list[float] = [0.0, 0.0, 1.0],
-        actor_residual_magnetic_field: list[float] = [0.0, 0.0, 0.0],
+        actor_residual_magnetic_field_body: list[float] = [0.0, 0.0, 0.0],
     ):
         """Add an attitude model to the actor based on initial conditions: attitude (roll, pitch & yaw angles)
         and angular velocity vector, modeling the evolution of the user specified pointing vector.
@@ -588,9 +588,8 @@ def set_attitude_model(
             actor (SpacecraftActor): Actor to model.
             actor_initial_attitude_in_rad (list of floats): Actor's initial attitude. Defaults to [0.0, 0.0, 0.0].
             actor_initial_angular_velocity (list of floats): Actor's initial angular velocity. Defaults to [0.0, 0.0, 0.0].
-            actor_pointing_vector_body (list of floats): Actor's pointing vector. Defaults to [0.0, 0.0, 1.0].
-            actor_residual_magnetic_field (list of floats): Actor's residual magnetic dipole moment vector.
-            Defaults to [0.0, 0.0, 0.0].
+            actor_pointing_vector_body (list of floats): Actor's pointing vector with respect to the body frame. Defaults to [0.0, 0.0, 1.0].
+            actor_residual_magnetic_field_body (list of floats): Actor's residual magnetic dipole moment vector in the body frame. Only needed if magnetic torque disturbances are modelled. Please, refer to [Tai L. Chow (2006) p. 148 - 149]. Defaults to [0.0, 0.0, 0.0].
         """
         # check for spacecraft actor
         assert isinstance(
@@ -609,7 +608,7 @@ def set_attitude_model(
             actor_initial_attitude_in_rad=actor_initial_attitude_in_rad,
             actor_initial_angular_velocity=actor_initial_angular_velocity,
             actor_pointing_vector_body=actor_pointing_vector_body,
-            actor_residual_magnetic_field=actor_residual_magnetic_field,
+            actor_residual_magnetic_field_body=actor_residual_magnetic_field_body,
         )
 
     @staticmethod

paseos/actors/spacecraft_actor.py

@@ -125,28 +125,6 @@ def attitude_disturbances(self) -> list:
         """
         return self._attitude_model._disturbances
 
-    @property
-    def body_moment_of_inertia(self) -> np.array:
-        """Gives the moment of inertia of the actor, assuming constant density.
-
-        Returns:
-            np.array: Mass moments of inertia for the actor
-
-        I is the moment of inertia, in the form of [[Ixx Ixy Ixz]
-                                                    [Iyx Iyy Iyx]
-                                                    [Izx Izy Izz]]
-        """
-        return self._spacecraft_body_model._body_moment_of_inertia
-
-    @property
-    def body_center_of_gravity(self) -> np.array:
-        """Gives the volumetric center of mass of the actor.
-
-        Returns:
-            np.array: Coordinates of the center of gravity of the mesh.
-        """
-        return self._spacecraft_body_model._body_center_of_gravity
-
     @property
     def temperature_in_K(self) -> float:
         """Returns the current temperature of the actor in K.
@@ -251,3 +229,25 @@ def angular_velocity(self):
             np.ndarray (owega_x, omega_y, omega_z).
         """
         return self._attitude_model._actor_angular_velocity_eci
+
+    @property
+    def body_moment_of_inertia_body(self) -> np.array:
+        """Gives the moment of inertia of the actor, assuming constant density, with respect to the body frame.
+
+        Returns:
+            np.array: Mass moments of inertia for the actor
+
+        I is the moment of inertia, in the form of [[Ixx Ixy Ixz]
+                                                    [Iyx Iyy Iyx]
+                                                    [Izx Izy Izz]]
+        """
+        return self._spacecraft_body_model._body_moment_of_inertia_body
+
+    @property
+    def body_center_of_gravity_body(self) -> np.array:
+        """Gives the volumetric center of mass of the actor with respect to the body frame.
+
+        Returns:
+            np.array: Coordinates of the center of gravity of the mesh.
+        """
+        return self._spacecraft_body_model._body_center_of_gravity_body

paseos/actors/spacecraft_body_model.py

@@ -9,8 +9,8 @@ class SpacecraftBodyModel:
     """
 
     _body_mesh = None
-    _body_center_of_gravity = None
-    _body_moment_of_inertia = None
+    _body_center_of_gravity_body = None
+    _body_moment_of_inertia_body = None
 
     def __init__(self, actor_mass, vertices=None, faces=None, scale=1) -> None:
         """Describes the geometry of the spacecraft and outputs relevant parameters related to the spacecraft body.
@@ -32,8 +32,8 @@ def __init__(self, actor_mass, vertices=None, faces=None, scale=1) -> None:
 
         self._body_mass = actor_mass
         self._body_mesh = self._create_body_mesh(vertices=vertices, faces=faces, scale=scale)
-        self._body_moment_of_inertia = self._body_mesh.moment_inertia * self._body_mass
-        self._body_center_of_gravity = self._body_mesh.center_mass
+        self._body_moment_of_inertia_body = self._body_mesh.moment_inertia * self._body_mass
+        self._body_center_of_gravity_body = self._body_mesh.center_mass
 
     def _create_body_mesh(self, vertices=None, faces=None, scale=1):
         """Creates the mesh of the satellite. If no vertices input is given, it defaults to a cuboid scaled by the
@@ -99,24 +99,4 @@ def body_mesh(self) -> np.array:
         """
         return self._body_mesh
 
-    @property
-    def body_moment_of_inertia(self) -> np.array:
-        """Gives the moment of inertia of the actor, assuming constant density.
-
-        Returns:
-            np.array: Mass moments of inertia for the actor
-
-        I is the moment of inertia, in the form of [[Ixx Ixy Ixz]
-                                                    [Iyx Iyy Iyx]
-                                                    [Izx Izy Izz]]
-        """
-        return self._body_moment_of_inertia
 
-    @property
-    def body_center_of_gravity(self) -> np.array:
-        """Gives the volumetric center of mass of the actor.
-
-        Returns:
-            np.array: Coordinates of the center of gravity of the mesh.
-        """
-        return self._body_center_of_gravity

paseos/attitude/attitude_model.py

@@ -1,5 +1,6 @@
 import numpy as np
 import pykep as pk
+from ..actors.spacecraft_actor import SpacecraftActor
 
 from .disturbance_torques_utils import (
     compute_aerodynamic_torque,
@@ -16,6 +17,15 @@
     rpy_to_body,
 )
 
+from enum import Enum
+
+class TorqueDisturbanceModel(Enum):
+    """Describes the different torque disturbance models supported.
+    1 - SpacePlot
+    """
+    Aerodynamic = 1
+    Gravitational = 2
+    Magnetic = 3
 
 class AttitudeModel:
     """This model describes the attitude (Roll, Pitch and Yaw) evolution of a spacecraft actor.
@@ -51,7 +61,7 @@ def __init__(
         actor_initial_angular_velocity: list[float] = [0.0, 0.0, 0.0],
         # pointing vector in body frame: (defaults to body z-axis)
         actor_pointing_vector_body: list[float] = [0.0, 0.0, 1.0],
-        actor_residual_magnetic_field: list[float] = [0.0, 0.0, 0.0],
+        actor_residual_magnetic_field_body: list[float] = [0.0, 0.0, 0.0],
     ):
         """Creates an attitude model to model actor attitude based on
         initial conditions (initial attitude and angular velocity) and
@@ -64,9 +74,28 @@ def __init__(
             actor_initial_angular_velocity (list of floats): Actor's initial angular velocity vector.
                 Defaults to [0., 0., 0.].
             actor_pointing_vector_body (list of floats): User defined vector in the Actor body. Defaults to [0., 0., 1]
-            actor_residual_magnetic_field (list of floats): Actor's own magnetic field modeled as dipole moment vector.
+            actor_residual_magnetic_field_body (list of floats): Actor's own magnetic field modeled as dipole moment vector.
                 Defaults to [0., 0., 0.].
         """
+        assert (isinstance(local_actor, SpacecraftActor)
+            ), "local_actor must be a ""SpacecraftActor""."
+
+        assert (
+                np.asarray(actor_initial_attitude_in_rad).shape[0] == 3 and actor_initial_attitude_in_rad.ndim == 1
+            ), "actor_initial_attitude_in_rad shall be [3] shaped."
+
+        assert (
+                np.asarray(actor_initial_angular_velocity).shape[0] == 3 and actor_initial_angular_velocity.ndim == 1
+            ), "actor_initial_angular_velocity shall be [3] shaped."
+
+        assert (
+                np.asarray(actor_pointing_vector_body).shape[0] == 3 and actor_pointing_vector_body.ndim == 1
+            ), "actor_pointing_vector_body shall be [3] shaped."
+
+        assert (
+                np.asarray(actor_residual_magnetic_field_body).shape[0] == 3 and actor_residual_magnetic_field_body.ndim == 1
+            ), "actor_residual_magnetic_field_body shall be [3] shaped."
+
         self._actor = local_actor
         # convert to np.ndarray
         self._actor_attitude_in_rad = np.array(actor_initial_attitude_in_rad)
@@ -92,7 +121,7 @@ def __init__(
         )
 
         # actor residual magnetic field (modeled as dipole)
-        self._actor_residual_magnetic_field = np.array(actor_residual_magnetic_field)
+        self._actor_residual_magnetic_field_body = np.array(actor_residual_magnetic_field_body)
 
     def _nadir_vector(self):
         """Compute unit vector pointing towards earth, in the inertial frame.
@@ -120,15 +149,15 @@ def compute_disturbance_torque(self, position, velocity, euler_angles, current_t
 
         if self._disturbances is not None:
             # TODO add solar disturbance
-            if "aerodynamic" in self._actor.attitude_disturbances:
+            if TorqueDisturbanceModel.Aerodynamic in self._actor.attitude_disturbances:
                 T += compute_aerodynamic_torque(
                     position,
                     velocity,
                     self._actor.geometric_model.mesh,
                     self.actor_attitude_in_rad,
                     current_temperature_K,
                 )
-            if "gravitational" in self._actor.attitude_disturbances:
+            if TorqueDisturbanceModel.Gravitational in self._actor.attitude_disturbances:
                 # Extract nadir vectors in different reference systems
                 nadir_vector_in_rpy = eci_to_rpy(self._nadir_vector(), position, velocity)
                 nadir_vector_in_body = rpy_to_body(nadir_vector_in_rpy, euler_angles)
@@ -142,14 +171,14 @@ def compute_disturbance_torque(self, position, velocity, euler_angles, current_t
                     self._actor.central_body.planet,
                     nadir_vector_in_body,
                     earth_rotation_vector_in_body,
-                    self._actor.body_moment_of_inertia,
+                    self._actor.body_moment_of_inertia_body,
                     np.linalg.norm(position),
                 )
-            if "magnetic" in self._actor.attitude_disturbances:
+            if TorqueDisturbanceModel.Magnetic in self._actor.attitude_disturbances:
                 time = self._actor.local_time
                 T += compute_magnetic_torque(
                     m_earth=self._actor.central_body.magnetic_dipole_moment(time),
-                    m_sat=self._actor_residual_magnetic_field,
+                    m_sat=self._actor_residual_magnetic_field_body,
                     position=self._actor.get_position(time),
                     velocity=self._actor.get_position_velocity(time)[1],
                     attitude=self._actor_attitude_in_rad,
@@ -165,7 +194,7 @@ def _calculate_angular_acceleration(self, current_temperature_K):
         # TODO in the future control torques could be added
         # Euler's equation for rigid body rotation: a = I^(-1) (T - w x (Iw))
         # with: a = angular acceleration, body_moment_of_inertia = inertia matrix, T = torque vector, w = angular velocity
-        self._actor_angular_acceleration = np.linalg.inv(self._actor.body_moment_of_inertia) @ (
+        self._actor_angular_acceleration = np.linalg.inv(self._actor.body_moment_of_inertia_body) @ (
             self.compute_disturbance_torque(
                 position=np.array(self._actor.get_position(self._actor.local_time)),
                 velocity=np.array(self._actor.get_position_velocity(self._actor.local_time)[1]),
@@ -174,7 +203,7 @@ def _calculate_angular_acceleration(self, current_temperature_K):
             )
             - np.cross(
                 self._actor_angular_velocity,
-                self._actor.body_moment_of_inertia @ self._actor_angular_velocity,
+                self._actor.body_moment_of_inertia_body @ self._actor_angular_velocity,
             )
         )
 
@@ -255,16 +284,13 @@ def update_attitude(self, dt, current_temperature_K):
         # time
         t = self._actor.local_time
 
-        # position
-        position = np.array(self._actor.get_position(t))
-
         # next position
         next_position = np.array(
             self._actor.get_position(pk.epoch(t.mjd2000 + dt * pk.SEC2DAY, "mjd2000"))
         )
 
-        # velocity
-        velocity = np.array(self._actor.get_position_velocity(self._actor.local_time)[1])
+        # position, velocity
+        position, velocity = np.array(self._actor.get_position_velocity(self._actor.local_time)[1])
 
         # Initial body vectors expressed in rpy frame: (x, y, z, custom pointing vector)
         xb_rpy, yb_rpy, zb_rpy, pointing_vector_rpy = self._body_axes_in_rpy()

paseos/attitude/disturbance_torques_utils.py

@@ -1,7 +1,5 @@
-"""this file contains functions that return attitude disturbance torque vectors expressed in the actor body frame"""
+"""this file contains functions that return attitude disturbance torque vectors expressed in the actor body frame."""
 
-# OUTPUT NUMPY ARRAYS
-# OUTPUT IN BODY FRAME
 import numpy as np
 from ..utils.reference_frame_transfer import (
     rpy_to_body,
@@ -10,7 +8,7 @@
 )
 
 
-def compute_aerodynamic_torque(r, v, mesh, actor_attitude_in_rad, current_spacecraft_temperature_K):
+def compute_aerodynamic_torque(position, velocity, mesh, actor_attitude_in_rad, current_spacecraft_temperature_K, central_body_radius_m):
     """Calculates the aerodynamic torque on the satellite.
     The model used is taken from "Roto-Translational Spacecraft Formation Control Using Aerodynamic Forces"; Ran. S,
     Jihe W., et al.; 2017. The mass density of the atmosphere is calculated from the best linear fit of the data
@@ -22,25 +20,26 @@ def compute_aerodynamic_torque(r, v, mesh, actor_attitude_in_rad, current_spacec
     temperature of the spacecraft and of the gas is low.
 
     Args:
-        r (np.array): distance from the satellite and the Earth's center of mass [km].
-        v (np.array): velocity of the satellite in ECI reference frame [m/s].
+        position (np.array): distance from the satellite and the Earth's center of mass [m].
+        velocity (np.array): velocity of the satellite in ECI reference frame [m/s].
         mesh (trimesh): mesh of the satellite from the geometric model.
         actor_attitude_in_rad (np.array): spacecraft actor in rad.
         current_spacecraft_temperature_K (float): current temperature in Kelvin.
+        central_body_radius_m (float): central body radius [m].
      Returns:
          T (np.array): torque vector in the spacecraft body frame.
     """
     # Constants for aerodynamic coefficients calculation
     temperature_gas = 1000  # [K]
     R = 8.314462  # Universal Gas Constant, [J/(K mol)]
-    altitude = np.linalg.norm(r) - 6371  # [km]
+    altitude = np.linalg.norm(position) - central_body_radius_m  # [km]
     density = 10 ** (
-        -(altitude + 1285) / 151
+        -(altitude + 1285e3) / 151e3
     )  # equation describing the best linear fit for the data, [kg/m^3]
-    molecular_speed_ratio_t = np.linalg.norm(v) / np.sqrt(
+    molecular_speed_ratio_t = np.linalg.norm(velocity) / np.sqrt(
         2 * R * temperature_gas
     )  # Used in the Cd and Cl calculation
-    molecular_speed_ratio_r = np.linalg.norm(v) / np.sqrt(
+    molecular_speed_ratio_r = np.linalg.norm(velocity) / np.sqrt(
         2 * R * current_spacecraft_temperature_K
     )  # Used in the Cd and Cl calculation
     accommodation_parameter = 0.85
@@ -53,7 +52,7 @@ def compute_aerodynamic_torque(r, v, mesh, actor_attitude_in_rad, current_spacec
     face_normals_rpy = np.dot(transformation_matrix_rpy_body, face_normals_sbf.T).T
 
     #  Get the velocity and transform it in the Roll Pitch Yaw frame. Get the unit vector associated with the latter
-    v_rpy = eci_to_rpy(v, r, v)
+    v_rpy = eci_to_rpy(velocity, position, velocity)
     unit_v_rpy = v_rpy / np.linalg.norm(v_rpy)
 
     #  Loop to get the normals, the areas, the angle of attack and the centroids of the faces with airflow, confronting them
@@ -160,15 +159,15 @@ def compute_aerodynamic_torque(r, v, mesh, actor_attitude_in_rad, current_spacec
 
         # Drag force on the plate [k]. Direction along the velocity vector.
         force_drag[k] = (
-            -0.5 * density * C_d * area_faces_airflow[k] * np.linalg.norm(v) ** 2 * unit_v_rpy
+            -0.5 * density * C_d * area_faces_airflow[k] * np.linalg.norm(velocity) ** 2 * unit_v_rpy
         )
         # Lift force on the plate [k]. Direction along the (v x n) x v direction, lift vector defined to be in that
         # direction. Intermediate step to get v x n.
         v_x_n_vector = np.cross(unit_v_rpy, normals_faces_with_airflow[k])
         not_norm_lift_vector = np.cross(v_x_n_vector, unit_v_rpy)
         lift_vector = not_norm_lift_vector / np.linalg.norm(not_norm_lift_vector)
         force_lift[k] = (
-            -0.5 * density * C_l * area_faces_airflow[k] * np.linalg.norm(v) ** 2 * lift_vector
+            -0.5 * density * C_l * area_faces_airflow[k] * np.linalg.norm(velocity) ** 2 * lift_vector
         )
 
         # Torque calculated as the product between the distance of the centroid from the geometric center of the

paseos/central_body/central_body.py

@@ -250,7 +250,7 @@ def magnetic_dipole_moment(
                 wgs84.latlon(pole_latitude_in_deg, pole_longitude_in_deg).at(t_skyfield).position.m
             )
             # Multiply geomagnetic pole position unit vector with dipole moment strength
-            magnetic_dipole_moment = (
+            magnetic_dipole_moment = -(
                 dipole_north_direction / np.linalg.norm(dipole_north_direction) * strength_in_Am2
             )
         else:

paseos/tests/attitude_disturbances_test.py

@@ -190,7 +190,7 @@ def flux_density_vector(actor, frame="eci"):
 
     # Now, align the body magnetic dipole with the local Earth magnetic flux density vector
     # Earth magnetic flux density vector at start position is approximately:
-    B0 = np.array([-3.18159529e-09, 1.02244882e-07, -3.72362170e-08])
+    B0 = np.array([3.18159529e-09, -1.02244882e-07, 3.72362170e-08])
 
     B_direction = B0 / np.linalg.norm(B0)