Merge pull request #202 from aidotse/attitude-magnetic-disturbance

Magnetic disturbances

paseos/actors/actor_builder.py

@@ -553,7 +553,8 @@ def set_attitude_model(
             actor: SpacecraftActor,
             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_pointing_vector_body: list[float] = [0.0, 0.0, 1.0],
+            actor_residual_magnetic_field: list[float] = [0.0, 0.0, 0.0],
 
     ):
         """Add an attitude model to the actor based on initial conditions: attitude (roll, pitch & yaw angles)
@@ -564,14 +565,16 @@ def set_attitude_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._attitude_model = AttitudeModel(
             local_actor=actor,
             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,
 
         )
 

paseos/actors/base_actor.py

@@ -355,7 +355,10 @@ def get_disturbances(self):
         Returns:
             list[string]: name of disturbances
         """
-        return self._disturbances
+        if self.has_attitude_disturbances:
+            return self._disturbances
+        else:
+            return "No disturbances"
 
     def is_in_line_of_sight(
         self,

paseos/attitude/attitude_model.py

@@ -1,6 +1,7 @@
 import numpy as np
 import pykep as pk
 
+
 from paseos.attitude.disturbance_calculations import (
     calculate_aero_torque,
     calculate_magnetic_torque,
@@ -43,9 +44,10 @@ def __init__(
         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.]
+        actor_pointing_vector_body: list[float] = [0., 0., 1.],
+        actor_residual_magnetic_field: list[float] = [0., 0., 0.],
     ):
-        """ Creates an attitude model to model actor attitude based on
+        """Creates an attitude model to model actor attitude based on
         initial conditions (initial attitude and angular velocity) and
         external disturbance torques.
 
@@ -56,6 +58,8 @@ 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.
+                Defaults to [0., 0., 0.].
         """
         self._actor = local_actor
         # convert to np.ndarray
@@ -81,6 +85,9 @@ def __init__(
             np.array(self._actor.get_position_velocity(self._actor.local_time)[1]),
         )
 
+        # actor residual magnetic field (modeled as dipole)
+        self._actor_residual_magnetic_field = np.array(actor_residual_magnetic_field)
+
     def _nadir_vector(self):
         """Compute unit vector pointing towards earth, in the inertial frame.
 
@@ -99,12 +106,20 @@ def _calculate_disturbance_torque(self):
         T = np.array([0.0, 0.0, 0.0])
 
         if self._actor.has_attitude_disturbances:
+            # TODO add solar disturbance
             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()
+                time = self._actor.local_time
+                T += calculate_magnetic_torque(
+                    m_earth=self._actor.central_body.magnetic_dipole_moment(time),
+                    m_sat=self._actor_residual_magnetic_field,
+                    position=self._actor.get_position(time),
+                    velocity=self._actor.get_position_velocity(time)[1],
+                    attitude=self._actor_attitude_in_rad,
+                )
         return T
 
     def _calculate_angular_acceleration(self):

paseos/attitude/disturbance_calculations.py

@@ -3,6 +3,7 @@
 # OUTPUT NUMPY ARRAYS
 # OUTPUT IN BODY FRAME
 import numpy as np
+from ..utils.reference_frame_transfer import rpy_to_body, eci_to_rpy
 
 
 def calculate_aero_torque():
@@ -19,8 +20,42 @@ def calculate_grav_torque():
     return np.array(T)
 
 
-def calculate_magnetic_torque():
-    # calculations for torques
-    # T must be in actor body fixed frame (to be discussed)
-    T = [0, 0, 0]
-    return np.array(T)
+def calculate_magnetic_torque(m_earth, m_sat, position, velocity, attitude):
+    """Calculates the external disturbance torque acting on the actor due to the magnetic field of the earth.
+    a dipole magnetic field flux density is described by the formula: B = μ0/(4πr³) * [3 r_hat(r_hat ⋅ m) − m]
+    With μ0 = 4 π 1e-7 H/m (vacuum permeability), r = actor distance from dipole center, r_hat = unit position vector,
+    and m the magnetic dipole moment vector of the Earth (magnitude in the year 2000 = 7.79 x 10²² Am²)
+
+    The disturbance torque is then calculated by: T = m_sat x B
+    With m_sat the (residual) magnetic dipole moment vector of the actor, magnitude usually between 0.1 - 20 Am² (SMAD)
+
+    https://en.wikipedia.org/wiki/Magnetic_dipole (or Chow, Tai L. (2006). Introduction to electromagnetic theory:
+    a modern perspective, used formular on p. 149)
+
+    Args:
+        m_earth (np.ndarray): magnetic dipole moment vector of the Earth in Am²
+        m_sat (np.ndarray): magnetic dipole moment vector of the actor, magnitude usually between 0.1 - 20 Am²
+        position (tuple or np.ndarray): actor position
+        velocity (tuple or np.ndarray): actor velocity (used for frame transformation)
+        attitude (np.ndarray): actor velocity (used for frame transformation)
+
+    Returns: Disturbance torque vector T (nd.array) in Nm in the actor body frame
+    """
+    # convert to np.ndarray
+    position = np.array(position)
+    velocity = np.array(velocity)
+
+    # actor distance
+    r = np.linalg.norm(position)
+    # actor unit position vector
+    r_hat = position / r
+
+    # magnetic field flux density at actor's position in Earth inertial frame
+    B = 1e-7 * (3 * np.dot(m_earth, r_hat) * r_hat - m_earth) / (r**3)
+
+    # transform field vector to body frame
+    B = rpy_to_body(eci_to_rpy(B, position, velocity), attitude)
+
+    # disturbance torque:
+    T = np.cross(m_sat, B)
+    return T

paseos/central_body/central_body.py

@@ -5,6 +5,8 @@
 from loguru import logger
 from pyquaternion import Quaternion
 from skspatial.objects import Sphere
+from skyfield.api import wgs84, load
+
 import pykep as pk
 import numpy as np
 
@@ -102,7 +104,9 @@ def blocks_sun(self, actor, t: pk.epoch, plot=False) -> bool:
 
         # Compute central body position in solar reference frame
         r_central_body_heliocentric, _ = np.array(self._planet.eph(t))
-        logger.trace("r_central_body_heliocentric is" + str(r_central_body_heliocentric))
+        logger.trace(
+            "r_central_body_heliocentric is" + str(r_central_body_heliocentric)
+        )
 
         # Compute satellite / actor position in solar reference frame
         r_sat_central_body_frame = np.array(actor.get_position(t))
@@ -126,7 +130,12 @@ def is_between_actors(self, actor_1, actor_2, t: pk.epoch, plot=False) -> bool:
         Returns:
             bool: True if the central body is between the two actors
         """
-        logger.debug("Computing line of sight between actors: " + str(actor_1) + " " + str(actor_2))
+        logger.debug(
+            "Computing line of sight between actors: "
+            + str(actor_1)
+            + " "
+            + str(actor_2)
+        )
         pos_1 = actor_1.get_position_velocity(t)
         pos_2 = actor_2.get_position_velocity(t)
 
@@ -153,7 +162,12 @@ def is_between_points(
         Returns:
             bool: True if the central body is between the two points
         """
-        logger.debug("Computing line of sight between points: " + str(point_1) + " " + str(point_2))
+        logger.debug(
+            "Computing line of sight between points: "
+            + str(point_1)
+            + " "
+            + str(point_2)
+        )
 
         point_1 = np.array(point_1)
         point_2 = np.array(point_2)
@@ -183,8 +197,12 @@ def is_between_points(
                 mesh_triangles=self._mesh[1],
             )
         else:
-            logger.error("No mesh or encompassing sphere provided. Cannot check visibility.")
-            raise ValueError("No mesh or encompassing sphere provided. Cannot check visibility.")
+            logger.error(
+                "No mesh or encompassing sphere provided. Cannot check visibility."
+            )
+            raise ValueError(
+                "No mesh or encompassing sphere provided. Cannot check visibility."
+            )
 
     def _apply_rotation(self, point, epoch: pk.epoch):
         """Applies the inverse rotation of the central body to the given point. This way
@@ -210,3 +228,57 @@ def _apply_rotation(self, point, epoch: pk.epoch):
         # Rotate point
         q = Quaternion(axis=self._rotation_axis, angle=angle)
         return q.rotate(point)
+
+    def magnetic_dipole_moment(
+        self,
+        epoch: pk.epoch,
+        strength_in_Am2=7.79e22,
+        pole_latitude_in_deg=79.6,
+        pole_longitude_in_deg=-71.6,
+    ):
+        """Returns the time-dependent magnetic dipole moment vector of central body.
+        Default values are for Earth. Earth dipole moment vector determined from the northern geomagnetic pole position
+        using skyfield api, and actor epoch. To model the simplified Earth magnetic field as a magnetic dipole with an
+        offset from the Earth rotational axis, at a specific point in time.
+
+        Earth geomagnetic pole position and dipole moment strength values from the year 2000:
+        Latitude: 79.6° N
+        Longitude: 71.6° W
+        Dipole moment: 7.79 x 10²² Am²
+        https://wdc.kugi.kyoto-u.ac.jp/poles/polesexp.html
+
+        (The same method used as ground station actor position determination)
+
+        Args:
+            epoch (pk.epoch): Epoch at which to get the dipole
+            strength_in_Am2 (float): dipole strength in Am². Defaults to 7.79e22
+            pole_latitude_in_deg (float): latitude of the Northern geomagnetic pole in degrees. Defaults to 79.6
+            pole_longitude_in_deg (float): longitude of the Northern geomagnetic pole in degrees. Defaults to -71.6
+
+        Returns: Time-dependent dipole moment vector in inertial frame (np.ndarray): [mx, my, mz]
+        """
+        if self.planet.name == "earth":
+            # Converting time to skyfield to use its API
+            t_skyfield = load.timescale().tt_jd(epoch.jd)
+
+            # North geomagnetic pole location on Earth surface in cartesian coordinates
+            dipole_north_direction = np.array(
+                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 = (
+                dipole_north_direction
+                / np.linalg.norm(dipole_north_direction)
+                * strength_in_Am2
+            )
+        else:
+            # TODO add other planets' magnetic fields
+            logger.warning(
+                "Magnetic dipole moment only modeled for Earth"
+            )
+            # For now: assume alignment with rotation axis
+            magnetic_dipole_moment = np.array([0, 0, 1]) * strength_in_Am2
+
+        return magnetic_dipole_moment