Attitude model v2

paseos/actors/actor_builder.py

@@ -532,7 +532,7 @@ def set_disturbances(
         Args:
             actor (SpacecraftActor): The actor to add to.
             aerodynamic (bool): Whether to consider aerodynamic disturbances in the attitude model. Defaults to False
-            gravitational (bool): Whether to consider gravitational disturbances in the attitude model. Defaults to False
+            gravitational (bool): Whether to consider gravity disturbances in the attitude model. Defaults to False
             magnetic (bool): Whether to consider magnetic disturbances in the attitude model. Defaults to False
         """
         # Create a list with user specified disturbances which are considered in the attitude modelling.
@@ -562,7 +562,7 @@ def set_attitude_model(
         Args:
             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_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].
         """
 

paseos/actors/base_actor.py

@@ -169,6 +169,15 @@ def has_attitude_model(self) -> bool:
         """
         return hasattr(self, "_attitude_model") and self._attitude_model is not None
 
+    @property
+    def has_attitude_disturbances(self) -> bool:
+        """Returns true if actor has attitude disturbances attributed, else false.
+
+        Returns:
+            bool: bool indicating presence.
+        """
+        return hasattr(self, "_disturbances") and self._disturbances is not None
+
     @property
     def mass(self) -> float:
         """Returns actor's mass in kg.

paseos/attitude/attitude_model.py

@@ -36,40 +36,50 @@ class AttitudeModel:
     _actor_pointing_vector_body = None
     _actor_pointing_vector_eci = None
 
-    # _actor_prev_pos = None
     def __init__(
         self,
         local_actor,
-        # initial conditions: (defaults to 0)
-        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],
-        # pointing vector in body frame: (defaults to z-axis)
-        actor_pointing_vector_body: list[float] = [0.0, 0.0, 1.0]
-        ## add args with default value = ...
-        # actor_dipole
-        # actor_drag_coefficient (more for geometric model)
-        # body_J2
+        # initial conditions:
+        actor_initial_attitude_in_rad: list[float] = [0., 0., 0.],
+        actor_initial_angular_velocity: list[float] = [0., 0., 0.],
+        # pointing vector in body frame: (defaults to body z-axis)
+        actor_pointing_vector_body: list[float] = [0., 0., 1.]
     ):
+        """ Creates an attitude model to model actor attitude based on
+        initial conditions (initial attitude and angular velocity) and
+        external disturbance torques.
+
+        Args:
+            actor (SpacecraftActor): Actor to model.
+            actor_initial_attitude_in_rad (list of floats): Actor's initial attitude ([roll, pitch, yaw]) angles.
+                Defaults to [0., 0., 0.].
+            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]
+        """
         self._actor = local_actor
+        # convert to np.ndarray
         self._actor_attitude_in_rad = np.array(actor_initial_attitude_in_rad)
         self._actor_angular_velocity = np.array(actor_initial_angular_velocity)
 
         # normalize inputted pointing vector & convert to np.ndarray
         self._actor_pointing_vector_body = np.array(
             actor_pointing_vector_body
         ) / np.linalg.norm(np.array(actor_pointing_vector_body))
+
         # pointing vector expressed in Earth-centered inertial frame
         self._actor_pointing_vector_eci = rpy_to_eci(
             body_to_rpy(self._actor_pointing_vector_body, self._actor_attitude_in_rad),
             np.array(self._actor.get_position(self._actor.local_time)),
             np.array(self._actor.get_position_velocity(self._actor.local_time)[1]),
         )
+
+        # angular velocity vector expressed in Earth-centered inertial frame
         self._actor_angular_velocity_eci = rpy_to_eci(
             body_to_rpy(self._actor_angular_velocity, self._actor_attitude_in_rad),
             np.array(self._actor.get_position(self._actor.local_time)),
             np.array(self._actor.get_position_velocity(self._actor.local_time)[1]),
         )
-        self._first_run = True
 
     def nadir_vector(self):
         # unused but might be useful during disturbance calculations or pointing vector relative position
@@ -87,22 +97,23 @@ def calculate_disturbance_torque(self):
         Returns:
             np.array([Tx, Ty, Tz]): total combined torques in Nm expressed in the spacecraft body frame.
         """
-
         T = np.array([0.0, 0.0, 0.0])
-        if "aerodynamic" in self._actor.get_disturbances():
-            T += calculate_aero_torque()
-        if "gravitational" in self._actor.get_disturbances():
-            T += calculate_grav_torque()
-        if "magnetic" in self._actor.get_disturbances():
-            T += calculate_magnetic_torque()
+
+        if self._actor.has_attitude_disturbances:
+            if "aerodynamic" in self._actor.get_disturbances():
+                T += calculate_aero_torque()
+            if "gravitational" in self._actor.get_disturbances():
+                T += calculate_grav_torque()
+            if "magnetic" in self._actor.get_disturbances():
+                T += calculate_magnetic_torque()
         return T
 
     def calculate_angular_acceleration(self):
         """Calculate the spacecraft angular acceleration (external disturbance torques and gyroscopic accelerations)."""
         # TODO in the future control torques could be added
 
         # moment of Inertia matrix:
-        I = self._actor._moment_of_inertia
+        I = self._actor._moment_of_inertia()
 
         # Euler's equation for rigid body rotation: a = I^(-1) (T - w x (Iw))
         # with: a = angular acceleration, I = inertia matrix, T = torque vector, w = angular velocity

paseos/tests/attitude_test.py

@@ -0,0 +1,170 @@
+"""Tests to see whether the attitude and disturbance models work as intended"""
+
+import numpy as np
+import pykep as pk
+import sys
+
+sys.path.append("../..")
+import paseos
+from paseos import ActorBuilder, SpacecraftActor
+
+
+def attitude_model_test():
+    """Testing the attitude model with no disturbances and known angular velocity.
+    One actor has orbit in Earth inertial x-y plane (equatorial) with initial velocity which rotates the actor with 180°
+    in 20 steps advancing time 100 seconds (pi/ (20 * 100)).
+    Another one has zero initial angular velocity.
+    This test mainly checks whether the model correctly models constant angular velocity without disturbances
+    """
+
+    earth = pk.planet.jpl_lp("earth")
+
+    # First actor constant angular acceleration
+    omega = np.pi/2000
+
+    # Define first local actor with angular velocity
+    sat1 = ActorBuilder.get_actor_scaffold("sat1", SpacecraftActor, pk.epoch(0))
+    ActorBuilder.set_orbit(sat1, [7000000, 0, 0], [0, 8000.0, 0], pk.epoch(0), earth)
+    ActorBuilder.set_geometric_model(sat1, mass=100)
+    ActorBuilder.set_attitude_model(
+        sat1,
+        actor_initial_angular_velocity=[0.0, omega, 0.0],
+        actor_pointing_vector_body=[0, 0, 1],
+    )
+
+    # Define second local actor without angular velocity
+    sat2 = ActorBuilder.get_actor_scaffold("sat2", SpacecraftActor, pk.epoch(0))
+    ActorBuilder.set_orbit(sat2, [7000000, 0, 0], [0, 8000.0, 0], pk.epoch(0), earth)
+    ActorBuilder.set_geometric_model(sat2, mass=100)
+    ActorBuilder.set_attitude_model(
+        sat2,
+        actor_initial_angular_velocity=[0.0, 0.0, 0.0],
+        actor_pointing_vector_body=[0, 0, 1],
+    )
+
+    # Check Initial values
+
+    # sat1
+    assert np.all(sat1._attitude_model._actor_pointing_vector_body == [0.0, 0.0, 1.0])
+    assert np.all(sat1._attitude_model._actor_pointing_vector_eci == [-1.0, 0.0, 0.0])
+    assert np.all(sat1._attitude_model._actor_attitude_in_rad == [0.0, 0.0, 0.0])
+    # positive angular velocity in body y direction is negative angular velocity in Earth inertial z direction:
+    assert np.all(sat1._attitude_model._actor_angular_velocity_eci == [0.0, 0.0, -omega])
+
+    # sat2
+    assert np.all(sat2._attitude_model._actor_pointing_vector_body == [0.0, 0.0, 1.0])
+    assert np.all(sat2._attitude_model._actor_pointing_vector_eci == [-1.0, 0.0, 0.0])
+    assert np.all(sat2._attitude_model._actor_attitude_in_rad == [0.0, 0.0, 0.0])
+
+    # Initialise simulation
+    sim = paseos.init_sim(sat1)
+    sim.add_known_actor(sat2)
+    # Run simulation 20 steps
+    for i in range(20):
+        sim.advance_time(100, 0)
+
+    # Testing the simulation went as intended
+    # Pointing vector from sat1 must be rotated from [-1, 0, 0] to [1, 0, 0]:
+    assert np.all(np.isclose(sat1.pointing_vector(), np.array([1.0, 0.0, 0.0])))
+    # Sat1 angular velocity in the body frame must stay constant:
+    assert np.all(
+        np.isclose(
+            sat1._attitude_model._actor_angular_velocity,
+            np.array([0.0, omega, 0.0]),
+        )
+    )
+    # Sat1 angular velocity in the Earth inertial frame must stay constant:
+    assert np.all(
+        np.isclose(
+            sat1.angular_velocity(),
+            np.array([0.0, 0.0, -omega]),
+        )
+    )
+
+    # Pointing vector from sat2 must not be rotated.
+    assert np.all(sat2.pointing_vector() == np.array([-1.0, 0.0, 0.0]))
+    # Sat2 angular velocity in the body frame must stay zero:
+    assert np.all(sat2._attitude_model._actor_angular_velocity == np.array([0.0, 0.0, 0.0]))
+
+
+def attitude_thermal_model_test():
+    """Testing the attitude model with no disturbances and no angular velocity, and ensuring the attitude model doesn't
+    break the thermal model (or vice versa)"""
+    earth = pk.planet.jpl_lp("earth")
+
+    # Define local actor
+    sat1 = ActorBuilder.get_actor_scaffold("sat1", SpacecraftActor, pk.epoch(0))
+    ActorBuilder.set_orbit(sat1, [7000000, 0, 0], [0, 8000.0, 0], pk.epoch(0), earth)
+    ActorBuilder.set_geometric_model(sat1, mass=100)
+    ActorBuilder.set_thermal_model(
+        actor=sat1,
+        actor_mass=sat1.mass,
+        actor_initial_temperature_in_K=273.15,
+        actor_sun_absorptance=1.0,
+        actor_infrared_absorptance=1.0,
+        actor_sun_facing_area=1.0,
+        actor_central_body_facing_area=1.0,
+        actor_emissive_area=1.0,
+        actor_thermal_capacity=1000,
+    )
+    ActorBuilder.set_attitude_model(sat1, actor_pointing_vector_body=[0, 0, 1])
+
+    # Check Initial values
+    assert np.all(sat1._attitude_model._actor_pointing_vector_body == [0.0, 0.0, 1.0])
+    assert np.all(sat1._attitude_model._actor_attitude_in_rad == [0.0, 0.0, 0.0])
+    assert sat1.temperature_in_K == 273.15
+
+    # Initialise simulation
+    sim = paseos.init_sim(sat1)
+
+    # Run simulation 20 steps
+    for i in range(21):
+        vector = sat1.pointing_vector()
+        sim.advance_time(100, 0)
+
+    # Testing the simulation went as intended
+    assert vector[0] == -1.0
+    assert sat1._attitude_model._actor_angular_velocity[1] == 0.0
+    assert np.round(sat1._attitude_model._actor_attitude_in_rad[0], 3) == 3.142
+    assert np.round(sat1.temperature_in_K, 3) == 278.522
+
+
+def attitude_and_orbit_test():
+    """This test checks both the orbit calculations, as well as the attitude.
+    The input is a simple orbit, and the angular velocity if 2pi/period. This means the initial conditions should be
+    the same as the conditions after one orbit"""
+
+    earth = pk.planet.jpl_lp("earth")
+
+    # Define local actor
+    sat1 = ActorBuilder.get_actor_scaffold("sat1", SpacecraftActor, pk.epoch(0))
+    ActorBuilder.set_orbit(sat1, [7000000, 0, 0], [0, 5460.0, 0], pk.epoch(0), earth)
+    ActorBuilder.set_geometric_model(sat1, mass=100)
+    orbit_period = 2 * np.pi * np.sqrt((6371000 + 7000000) ** 3 / 3.986004418e14)
+    ActorBuilder.set_attitude_model(
+        sat1,
+        actor_initial_angular_velocity=[0.0, 2 * np.pi / orbit_period, 0.0],
+        actor_pointing_vector_body=[0, 0, 1],
+    )
+    ActorBuilder.set_disturbances(sat1)
+
+    # Check Initial values
+    assert np.all(sat1._attitude_model._actor_pointing_vector_body == [0.0, 0.0, 1.0])
+    assert np.all(sat1._attitude_model._actor_attitude_in_rad == [0.0, 0.0, 0.0])
+    vector = sat1.pointing_vector()
+    assert vector[0] == -1.0
+
+    # Initialise simulation
+    sim = paseos.init_sim(sat1)
+
+    # Run simulation 10 steps
+    for i in range(11):
+        vector = sat1.pointing_vector()
+        sim.advance_time(orbit_period / 10, 0)
+
+    # Testing the simulation went as intended
+    assert sat1._attitude_model._actor_pointing_vector_body[2] == 1.0
+    assert vector[0] == -1.0
+
+
+attitude_model_test()