From 2f9969f04d6e370494937071091b7c7395cc51cf Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Thu, 25 Jan 2024 23:38:44 +0100
Subject: [PATCH 01/16] test commit for seeing if merge shows up.

---
 paseos/attitude/attitude_model.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index 1b278f7..b699550 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -95,6 +95,7 @@ def calculate_disturbance_torque(self):
             if "gravitational" in self._actor.get_disturbances():
                 T += calculate_grav_torque()
             if "magnetic" in self._actor.get_disturbances():
+                # add something like this: self.central_body.rotate_body()
                 T += calculate_magnetic_torque()
         return T
 

From b56c0765a1fecdddf45e6de433c838da87f987b6 Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Sat, 27 Jan 2024 13:17:54 +0100
Subject: [PATCH 02/16] added Earth magnetic dipole moment vector function

---
 paseos/actors/base_actor.py                 |  5 ++-
 paseos/attitude/attitude_model.py           | 39 +++++++++++++++++++--
 paseos/attitude/disturbance_calculations.py |  2 +-
 3 files changed, 42 insertions(+), 4 deletions(-)

diff --git a/paseos/actors/base_actor.py b/paseos/actors/base_actor.py
index 150cb5d..a5dc9e0 100644
--- a/paseos/actors/base_actor.py
+++ b/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,
diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index b699550..ccaec82 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -1,5 +1,7 @@
 import numpy as np
 import pykep as pk
+from skyfield.api import wgs84, load
+
 
 from paseos.attitude.disturbance_calculations import (
     calculate_aero_torque,
@@ -81,6 +83,38 @@ def nadir_vector(self):
         u = np.array(self._actor.get_position(self._actor.local_time))
         return -u / np.linalg.norm(u)
 
+    def earth_magnetic_dipole_moment(self):
+        """Returns the Earth magnetic dipole moment vector 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.
+
+        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²²
+        https://wdc.kugi.kyoto-u.ac.jp/poles/polesexp.html
+
+        (The same method used as ground station actor position determination)
+
+        Returns: Time-dependent Earth dipole moment vector in ECI (np.ndarray): [mx, my, mz]
+        """
+        # North geomagnetic pole location
+        lat = 79.6
+        lon = -71.6
+
+        # Converting time to skyfield to use its API
+        t_skyfield = load.timescale().tt_jd(self._actor.local_time.jd)
+
+        # North geomagnetic pole location on Earth surface in cartesian coordinates
+        dipole_north_direction = np.array(
+            wgs84.latlon(lat, lon).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) * 7.79e22
+        )
+        return magnetic_dipole_moment
+
     def calculate_disturbance_torque(self):
         """Compute total torque due to user specified disturbances.
 
@@ -95,8 +129,9 @@ def calculate_disturbance_torque(self):
             if "gravitational" in self._actor.get_disturbances():
                 T += calculate_grav_torque()
             if "magnetic" in self._actor.get_disturbances():
-                # add something like this: self.central_body.rotate_body()
-                T += calculate_magnetic_torque()
+                # Earth magnetic dipole moment at right position
+                magnetic_dipole_moment = self.earth_magnetic_dipole_moment()
+                T += calculate_magnetic_torque(magnetic_dipole_moment)
         return T
 
     def calculate_angular_acceleration(self):
diff --git a/paseos/attitude/disturbance_calculations.py b/paseos/attitude/disturbance_calculations.py
index 43f6fe2..f2ebd34 100644
--- a/paseos/attitude/disturbance_calculations.py
+++ b/paseos/attitude/disturbance_calculations.py
@@ -19,7 +19,7 @@ def calculate_grav_torque():
     return np.array(T)
 
 
-def calculate_magnetic_torque():
+def calculate_magnetic_torque(earth_m, sat_m, pos):
     # calculations for torques
     # T must be in actor body fixed frame (to be discussed)
     T = [0, 0, 0]

From ce745722f7ed1d6e27037f5d250fe33341c0643e Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Sat, 27 Jan 2024 13:38:07 +0100
Subject: [PATCH 03/16] Commit (unimportant or placeholder files/changes)

---
 paseos/attitude/attitude_model.py            |  4 +++-
 test_attitude_code/test_attitude_plotting.py | 18 ++++--------------
 2 files changed, 7 insertions(+), 15 deletions(-)

diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index ccaec82..dde97c1 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -131,7 +131,9 @@ def calculate_disturbance_torque(self):
             if "magnetic" in self._actor.get_disturbances():
                 # Earth magnetic dipole moment at right position
                 magnetic_dipole_moment = self.earth_magnetic_dipole_moment()
-                T += calculate_magnetic_torque(magnetic_dipole_moment)
+                m_sat = "placeholder"
+                position = "placeholder"
+                T += calculate_magnetic_torque(magnetic_dipole_moment, m_sat, position)
         return T
 
     def calculate_angular_acceleration(self):
diff --git a/test_attitude_code/test_attitude_plotting.py b/test_attitude_code/test_attitude_plotting.py
index 46ed59e..09a56e2 100644
--- a/test_attitude_code/test_attitude_plotting.py
+++ b/test_attitude_code/test_attitude_plotting.py
@@ -25,18 +25,6 @@
 )
 ActorBuilder.set_geometric_model(sat1, mass=500)
 
-ActorBuilder.set_thermal_model(
-    sat1,
-    actor_mass=100,
-    actor_initial_temperature_in_K=270,
-    actor_sun_absorptance=0.5,
-    actor_infrared_absorptance=0.5,
-    actor_sun_facing_area=1,
-    actor_central_body_facing_area=4,
-    actor_emissive_area=18,
-    actor_thermal_capacity=0.89,
-)
-
 # when i = 21 in loop and advance time =100, pi/2000 rad/sec will rotate 180 deg about 1 axis
 ActorBuilder.set_attitude_model(
     sat1,
@@ -45,7 +33,7 @@
     actor_initial_attitude_in_rad=[0.0, 0.0, 0.0],
 )
 # disturbances:
-# ActorBuilder.set_disturbances(sat1,True, True)
+ActorBuilder.set_disturbances(sat1, False, False, True)
 # ActorBuilder.set_disturbances(sat1)
 
 sim = paseos.init_sim(sat1)
@@ -88,11 +76,12 @@
     print(
         "plotted attitude:", euler, " at position: ", pos, " pointing v: ", vector / 2e6
     )
+    m = sat1._attitude_model.earth_magnetic_dipole_moment() * 1e-16
     # plot vectors
     ax.quiver(pos[0], pos[1], pos[2], ang_vel[0], ang_vel[1], ang_vel[2], color="m")
     ax.quiver(pos[0], pos[1], pos[2], vector[0], vector[1], vector[2])
+    ax.quiver(0, 0, 0, m[0], m[1], m[2], color="g")
 
-    # sim.advance_time(100, 0)
     sim.advance_time(100, 0)
 
 
@@ -111,3 +100,4 @@
 ax.plot(x, y, z)
 ax.scatter(0, 0, 0)
 plt.show()
+

From bdcd11b1ab89446853eb69676c6342e3d9da3b28 Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Sat, 27 Jan 2024 21:27:54 +0100
Subject: [PATCH 04/16] magnetic disturbance code v1

---
 paseos/attitude/attitude_model.py           | 24 ++++++------
 paseos/attitude/disturbance_calculations.py | 42 ++++++++++++++++++---
 2 files changed, 47 insertions(+), 19 deletions(-)

diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index dde97c1..b3d5c6a 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -38,7 +38,7 @@ class AttitudeModel:
     _actor_pointing_vector_body = None
     _actor_pointing_vector_eci = None
 
-    # _actor_prev_pos = None
+
     def __init__(
         self,
         local_actor,
@@ -46,11 +46,8 @@ def __init__(
         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
+        actor_pointing_vector_body: list[float] = [0.0, 0.0, 1.0],
+        actor_residual_magnetic_field: list[float] = None
     ):
         self._actor = local_actor
         self._actor_attitude_in_rad = np.array(actor_initial_attitude_in_rad)
@@ -71,7 +68,7 @@ def __init__(
             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
+        self._actor_residual_magnetic_field = np.array(actor_residual_magnetic_field)
 
     def nadir_vector(self):
         # unused but might be useful during disturbance calculations or pointing vector relative position
@@ -91,7 +88,7 @@ def earth_magnetic_dipole_moment(self):
         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²²
+        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)
@@ -129,11 +126,12 @@ def calculate_disturbance_torque(self):
             if "gravitational" in self._actor.get_disturbances():
                 T += calculate_grav_torque()
             if "magnetic" in self._actor.get_disturbances():
-                # Earth magnetic dipole moment at right position
-                magnetic_dipole_moment = self.earth_magnetic_dipole_moment()
-                m_sat = "placeholder"
-                position = "placeholder"
-                T += calculate_magnetic_torque(magnetic_dipole_moment, m_sat, position)
+                T += calculate_magnetic_torque(
+                    m_earth=self.earth_magnetic_dipole_moment(),
+                    m_sat=self._actor_residual_magnetic_field,
+                    position=self._actor.get_position(self._actor.local_time),
+                    velocity=self._actor.get_position_velocity(self._actor.local_time)[1],
+                    attitude=self._actor_attitude_in_rad)
         return T
 
     def calculate_angular_acceleration(self):
diff --git a/paseos/attitude/disturbance_calculations.py b/paseos/attitude/disturbance_calculations.py
index f2ebd34..f9b9073 100644
--- a/paseos/attitude/disturbance_calculations.py
+++ b/paseos/attitude/disturbance_calculations.py
@@ -3,7 +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():
     # calculations for torques
@@ -19,8 +19,38 @@ def calculate_grav_torque():
     return np.array(T)
 
 
-def calculate_magnetic_torque(earth_m, sat_m, pos):
-    # 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 (np.ndarray): actor position
+        velocity (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
+    """
+    # actor distance and unit position vector
+    r = np.linalg.norm(position)
+    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

From f64cc0eb82f70f2c21a00e6617fad520b7742e5b Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Sat, 27 Jan 2024 22:43:58 +0100
Subject: [PATCH 05/16] fixed some stuff, added way to visualize B on actor
 position (needs to be deleted)

---
 paseos/actors/actor_builder.py               |  7 ++++--
 paseos/attitude/attitude_model.py            | 24 +++++++++++++++-----
 paseos/attitude/disturbance_calculations.py  | 14 ++++++++----
 test_attitude_code/test_attitude_plotting.py | 15 ++++++++----
 4 files changed, 42 insertions(+), 18 deletions(-)

diff --git a/paseos/actors/actor_builder.py b/paseos/actors/actor_builder.py
index c567ecf..f6b0d8c 100644
--- a/paseos/actors/actor_builder.py
+++ b/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,6 +565,8 @@ 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(
@@ -571,7 +574,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,
 
         )
 
diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index b3d5c6a..46ade23 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -38,7 +38,6 @@ class AttitudeModel:
     _actor_pointing_vector_body = None
     _actor_pointing_vector_eci = None
 
-
     def __init__(
         self,
         local_actor,
@@ -47,7 +46,7 @@ def __init__(
         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],
-        actor_residual_magnetic_field: list[float] = None
+        actor_residual_magnetic_field: list[float] = [0.0, 0.0, 0.0],
     ):
         self._actor = local_actor
         self._actor_attitude_in_rad = np.array(actor_initial_attitude_in_rad)
@@ -69,6 +68,7 @@ def __init__(
             np.array(self._actor.get_position_velocity(self._actor.local_time)[1]),
         )
         self._actor_residual_magnetic_field = np.array(actor_residual_magnetic_field)
+        self._actor_magnetic_flux = np.array([0, 0, 0])
 
     def nadir_vector(self):
         # unused but might be useful during disturbance calculations or pointing vector relative position
@@ -119,6 +119,12 @@ def calculate_disturbance_torque(self):
             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])
+        dt = 10
+        # time
+        t = self._actor.local_time
+        next_position = np.array(
+            self._actor.get_position(pk.epoch(t.mjd2000 + dt * pk.SEC2DAY, "mjd2000"))
+        )
 
         if self._actor.has_attitude_disturbances:
             if "aerodynamic" in self._actor.get_disturbances():
@@ -126,12 +132,18 @@ def calculate_disturbance_torque(self):
             if "gravitational" in self._actor.get_disturbances():
                 T += calculate_grav_torque()
             if "magnetic" in self._actor.get_disturbances():
-                T += calculate_magnetic_torque(
+                t, b = calculate_magnetic_torque(
                     m_earth=self.earth_magnetic_dipole_moment(),
                     m_sat=self._actor_residual_magnetic_field,
-                    position=self._actor.get_position(self._actor.local_time),
-                    velocity=self._actor.get_position_velocity(self._actor.local_time)[1],
-                    attitude=self._actor_attitude_in_rad)
+                    # position=self._actor.get_position(self._actor.local_time),
+                    position=next_position,
+                    velocity=self._actor.get_position_velocity(self._actor.local_time)[
+                        1
+                    ],
+                    attitude=self._actor_attitude_in_rad,
+                )
+                T += t
+                self._actor_magnetic_flux = b
         return T
 
     def calculate_angular_acceleration(self):
diff --git a/paseos/attitude/disturbance_calculations.py b/paseos/attitude/disturbance_calculations.py
index f9b9073..4726d44 100644
--- a/paseos/attitude/disturbance_calculations.py
+++ b/paseos/attitude/disturbance_calculations.py
@@ -34,23 +34,27 @@ def calculate_magnetic_torque(m_earth, m_sat, position, velocity, attitude):
     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 (np.ndarray): actor position
-        velocity (np.ndarray): actor velocity (used for frame transformation)
+        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 and unit position vector
     r = np.linalg.norm(position)
     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)
-
+    B_test = B
     # 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
+    print("B ", B, " x m ", m_sat, " = ", T)
+    return T, B_test
diff --git a/test_attitude_code/test_attitude_plotting.py b/test_attitude_code/test_attitude_plotting.py
index 09a56e2..90bdcf7 100644
--- a/test_attitude_code/test_attitude_plotting.py
+++ b/test_attitude_code/test_attitude_plotting.py
@@ -15,11 +15,14 @@
 
 # Define spacecraft actor
 sat1 = ActorBuilder.get_actor_scaffold("sat1", SpacecraftActor, pk.epoch(0))
+R = 3000000+6371000
+theta = 11*np.pi/180
 
 ActorBuilder.set_orbit(
     sat1,
-    position=[10000000, 0, 0],
-    velocity=[0, 8000, 0],
+    position=[R * np.cos(theta), 0, -R * np.sin(theta)],
+    #velocity=[0, 8000, 0],
+    velocity=[0, 6519.49, 0],
     epoch=pk.epoch(0),
     central_body=earth,
 )
@@ -28,9 +31,10 @@
 # when i = 21 in loop and advance time =100, pi/2000 rad/sec will rotate 180 deg about 1 axis
 ActorBuilder.set_attitude_model(
     sat1,
-    actor_initial_angular_velocity=[0.0, np.pi / 2000, 0.0],
+    actor_initial_angular_velocity=[0.0, 0 * np.pi / 2000, 0.0],
     actor_pointing_vector_body=[0.0, 0.0, 1.0],
     actor_initial_attitude_in_rad=[0.0, 0.0, 0.0],
+    actor_residual_magnetic_field=[0.0, 0.0, 0.05],
 )
 # disturbances:
 ActorBuilder.set_disturbances(sat1, False, False, True)
@@ -77,12 +81,14 @@
         "plotted attitude:", euler, " at position: ", pos, " pointing v: ", vector / 2e6
     )
     m = sat1._attitude_model.earth_magnetic_dipole_moment() * 1e-16
+    B = sat1._attitude_model._actor_magnetic_flux * 1e12
     # plot vectors
     ax.quiver(pos[0], pos[1], pos[2], ang_vel[0], ang_vel[1], ang_vel[2], color="m")
     ax.quiver(pos[0], pos[1], pos[2], vector[0], vector[1], vector[2])
     ax.quiver(0, 0, 0, m[0], m[1], m[2], color="g")
+    ax.quiver(pos[0], pos[1], pos[2], B[0], B[1], B[2], color="y")
 
-    sim.advance_time(100, 0)
+    sim.advance_time(300, 0)
 
 
 # 3D figure limits
@@ -100,4 +106,3 @@
 ax.plot(x, y, z)
 ax.scatter(0, 0, 0)
 plt.show()
-

From 30d108f80ffa72a9fa0ceb0ac04f2078fd9d0f51 Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Sun, 28 Jan 2024 18:46:42 +0100
Subject: [PATCH 06/16] next position in disturbance calculations and cleanup

---
 paseos/attitude/attitude_model.py           | 13 ++++---------
 paseos/attitude/disturbance_calculations.py | 11 ++++++-----
 2 files changed, 10 insertions(+), 14 deletions(-)

diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index 46ade23..abf9fdf 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -68,7 +68,6 @@ def __init__(
             np.array(self._actor.get_position_velocity(self._actor.local_time)[1]),
         )
         self._actor_residual_magnetic_field = np.array(actor_residual_magnetic_field)
-        self._actor_magnetic_flux = np.array([0, 0, 0])
 
     def nadir_vector(self):
         # unused but might be useful during disturbance calculations or pointing vector relative position
@@ -125,6 +124,7 @@ def calculate_disturbance_torque(self):
         next_position = np.array(
             self._actor.get_position(pk.epoch(t.mjd2000 + dt * pk.SEC2DAY, "mjd2000"))
         )
+        position = np.array(self._actor.get_position(t))
 
         if self._actor.has_attitude_disturbances:
             if "aerodynamic" in self._actor.get_disturbances():
@@ -132,18 +132,13 @@ def calculate_disturbance_torque(self):
             if "gravitational" in self._actor.get_disturbances():
                 T += calculate_grav_torque()
             if "magnetic" in self._actor.get_disturbances():
-                t, b = calculate_magnetic_torque(
+                T += calculate_magnetic_torque(
                     m_earth=self.earth_magnetic_dipole_moment(),
                     m_sat=self._actor_residual_magnetic_field,
-                    # position=self._actor.get_position(self._actor.local_time),
                     position=next_position,
-                    velocity=self._actor.get_position_velocity(self._actor.local_time)[
-                        1
-                    ],
+                    velocity=self._actor.get_position_velocity(self._actor.local_time)[1],
                     attitude=self._actor_attitude_in_rad,
                 )
-                T += t
-                self._actor_magnetic_flux = b
         return T
 
     def calculate_angular_acceleration(self):
@@ -151,7 +146,7 @@ def calculate_angular_acceleration(self):
         # 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
diff --git a/paseos/attitude/disturbance_calculations.py b/paseos/attitude/disturbance_calculations.py
index 4726d44..7e52795 100644
--- a/paseos/attitude/disturbance_calculations.py
+++ b/paseos/attitude/disturbance_calculations.py
@@ -5,6 +5,7 @@
 import numpy as np
 from ..utils.reference_frame_transfer import rpy_to_body, eci_to_rpy
 
+
 def calculate_aero_torque():
     # calculations for torques
     # T must be in actor body fixed frame (to be discussed)
@@ -44,17 +45,17 @@ def calculate_magnetic_torque(m_earth, m_sat, position, velocity, attitude):
     position = np.array(position)
     velocity = np.array(velocity)
 
-    # actor distance and unit position vector
+    # 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)
-    B_test = B
+    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)
-    print("B ", B, " x m ", m_sat, " = ", T)
-    return T, B_test
+    return T

From e95504bb484a993f0972d4aa54fb359394aa0a89 Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Sun, 28 Jan 2024 18:49:37 +0100
Subject: [PATCH 07/16] commit to save, but won't work (intentional)

---
 test_attitude_code/test_attitude_plotting.py | 83 ++++++++++++++------
 1 file changed, 60 insertions(+), 23 deletions(-)

diff --git a/test_attitude_code/test_attitude_plotting.py b/test_attitude_code/test_attitude_plotting.py
index 90bdcf7..4168231 100644
--- a/test_attitude_code/test_attitude_plotting.py
+++ b/test_attitude_code/test_attitude_plotting.py
@@ -15,32 +15,60 @@
 
 # Define spacecraft actor
 sat1 = ActorBuilder.get_actor_scaffold("sat1", SpacecraftActor, pk.epoch(0))
-R = 3000000+6371000
-theta = 11*np.pi/180
+sat2 = ActorBuilder.get_actor_scaffold("sat2", SpacecraftActor, pk.epoch(0))
+# geostationary orbit:
+
+lat = 79.6 * np.pi / 180
+lon = -71.6 * np.pi / 180
+
+R = 6371000 + 35786000
+v = 3074.66
+
+arr = np.array([-1.24217547e-09, 1.01735657e-07, -3.87201400e-08])
+initial_magn = np.ndarray.tolist(arr / np.linalg.norm(arr) * 20)
+initial_pointing = np.ndarray.tolist(arr / np.linalg.norm(arr))
 
 ActorBuilder.set_orbit(
     sat1,
-    position=[R * np.cos(theta), 0, -R * np.sin(theta)],
-    #velocity=[0, 8000, 0],
-    velocity=[0, 6519.49, 0],
+    position=[R * np.cos(np.pi / 2 + lon), R * np.sin(np.pi / 2 + lon), 0],
+    velocity=[-v * np.sin(np.pi / 2 + lon), v * np.cos(np.pi / 2 + lon), 0],
+    epoch=pk.epoch(0),
+    central_body=earth,
+)
+ActorBuilder.set_orbit(
+    sat2,
+    position=[R * np.cos(np.pi / 2 + lon), R * np.sin(np.pi / 2 + lon), 0],
+    velocity=[-v * np.sin(np.pi / 2 + lon), v * np.cos(np.pi / 2 + lon), 0],
     epoch=pk.epoch(0),
     central_body=earth,
 )
-ActorBuilder.set_geometric_model(sat1, mass=500)
+
+ActorBuilder.set_geometric_model(sat1, mass=100)
+ActorBuilder.set_geometric_model(sat2, mass=100)
 
 # when i = 21 in loop and advance time =100, pi/2000 rad/sec will rotate 180 deg about 1 axis
 ActorBuilder.set_attitude_model(
     sat1,
-    actor_initial_angular_velocity=[0.0, 0 * np.pi / 2000, 0.0],
-    actor_pointing_vector_body=[0.0, 0.0, 1.0],
+    actor_initial_angular_velocity=[0.0, 0.0, 0.0],
+    actor_pointing_vector_body=initial_pointing,
     actor_initial_attitude_in_rad=[0.0, 0.0, 0.0],
-    actor_residual_magnetic_field=[0.0, 0.0, 0.05],
+    actor_residual_magnetic_field=initial_magn,
 )
+ActorBuilder.set_attitude_model(
+    sat2,
+    actor_initial_angular_velocity=[0.0, 0.0, 0.0],
+    actor_pointing_vector_body=initial_pointing,
+    actor_initial_attitude_in_rad=[0.0, 0.0, 0.0],
+    actor_residual_magnetic_field=[0.0, 0.0, 0.0],
+)
+
 # disturbances:
-ActorBuilder.set_disturbances(sat1, False, False, True)
-# ActorBuilder.set_disturbances(sat1)
+ActorBuilder.set_disturbances(sat1, magnetic=True)
+ActorBuilder.set_disturbances(sat2, magnetic=True)
+
 
 sim = paseos.init_sim(sat1)
+sim.add_known_actor(sat2)
 plt.close()
 
 pos = []
@@ -51,7 +79,7 @@
 
 fig = plt.figure()
 ax = fig.add_subplot(111, projection="3d")
-for i in range(21):
+for i in range(46):
     print("----------", i, "----------")
     pos = sat1.get_position(sat1.local_time)
     x.append(sat1.get_position(sat1.local_time)[0])
@@ -64,7 +92,13 @@
     vector = sat1.pointing_vector()
     vector[np.isclose(vector, np.zeros(3))] = 0
     # scale for plotting
-    vector = vector * 2e6
+    vector = vector * 8e6
+
+    # pointing vector:
+    vector2 = sat2.pointing_vector()
+    vector2[np.isclose(vector2, np.zeros(3))] = 0
+    # scale for plotting
+    vector2 = vector2 * 8e6
 
     # angular velocity vector:
     # normalize first:
@@ -75,20 +109,23 @@
         ang_vel = sat1.angular_velocity() / np.linalg.norm(sat1.angular_velocity())
         ang_vel[np.isclose(ang_vel, np.zeros(3))] = 0
     # scale for plotting
-    ang_vel = ang_vel * 2e6
+    ang_vel = ang_vel * 2e7
+
+    # print("plotted attitude:", euler, " at position: ", pos, " pointing v: ", vector / 2e6)
 
-    print(
-        "plotted attitude:", euler, " at position: ", pos, " pointing v: ", vector / 2e6
-    )
-    m = sat1._attitude_model.earth_magnetic_dipole_moment() * 1e-16
-    B = sat1._attitude_model._actor_magnetic_flux * 1e12
+    m = sat1._attitude_model.earth_magnetic_dipole_moment() * 6e-16
+    B = sat1._attitude_model._actor_magnetic_flux * 2e14
+    print(sat1._attitude_model._actor_magnetic_flux)
     # plot vectors
-    ax.quiver(pos[0], pos[1], pos[2], ang_vel[0], ang_vel[1], ang_vel[2], color="m")
-    ax.quiver(pos[0], pos[1], pos[2], vector[0], vector[1], vector[2])
-    ax.quiver(0, 0, 0, m[0], m[1], m[2], color="g")
+    # ax.quiver(pos[0], pos[1], pos[2], ang_vel[0], ang_vel[1], ang_vel[2], color="m")
+    ax.quiver(pos[0], pos[1], pos[2], vector[0], vector[1], vector[2], linewidths=3, color="r")
+    ax.quiver(pos[0], pos[1], pos[2], vector2[0], vector2[1], vector2[2], linewidths=3, color="m")
+    if not i % 10:
+        ax.quiver(0, 0, 0, m[0], m[1], m[2], color="y")
     ax.quiver(pos[0], pos[1], pos[2], B[0], B[1], B[2], color="y")
 
-    sim.advance_time(300, 0)
+    sim.advance_time(1000, 0)
+    # sim.advance_time(3446, 0)
 
 
 # 3D figure limits

From 437b8d914616be49ca177ca1605331738bf77fcd Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Mon, 29 Jan 2024 15:13:06 +0100
Subject: [PATCH 08/16] test for magnetic disturbance

---
 paseos/tests/attitude_test.py                | 141 ++++++++++++++++++-
 test_attitude_code/test_attitude_plotting.py |  23 ++-
 2 files changed, 153 insertions(+), 11 deletions(-)

diff --git a/paseos/tests/attitude_test.py b/paseos/tests/attitude_test.py
index 25ef4ca..37acb28 100644
--- a/paseos/tests/attitude_test.py
+++ b/paseos/tests/attitude_test.py
@@ -7,6 +7,7 @@
 sys.path.append("../..")
 import paseos
 from paseos import ActorBuilder, SpacecraftActor
+from paseos.utils.reference_frame_transfer import eci_to_rpy, rpy_to_body
 
 
 def attitude_model_test():
@@ -20,7 +21,7 @@ def attitude_model_test():
     earth = pk.planet.jpl_lp("earth")
 
     # First actor constant angular acceleration
-    omega = np.pi/2000
+    omega = np.pi / 2000
 
     # Define first local actor with angular velocity
     sat1 = ActorBuilder.get_actor_scaffold("sat1", SpacecraftActor, pk.epoch(0))
@@ -49,7 +50,9 @@ def attitude_model_test():
     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])
+    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])
@@ -84,7 +87,9 @@ def attitude_model_test():
     # 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]))
+    assert np.all(
+        sat2._attitude_model._actor_angular_velocity == np.array([0.0, 0.0, 0.0])
+    )
 
 
 def attitude_thermal_model_test():
@@ -167,4 +172,132 @@ def attitude_and_orbit_test():
     assert vector[0] == -1.0
 
 
-attitude_model_test()
+def magnetic_disturbance_test():
+    """Tests the magnetic disturbance torques applied in the attitude model.
+    First: put two spacecraft actors in a geostationary orbit (disregarding the relative magnetic field rotation of the
+    Earth). Both actor's own magnetic dipole moment aligned with the local magnetic flux density vector of the Earth
+    magnetic field. One is non-magnetic and is expected to have a fixed attitude in the Earth inertial frame.
+    The other (magnetic) actor should stay aligned with the Earth magnetic field.
+    """
+    # Define central body
+    earth = pk.planet.jpl_lp("earth")
+
+    # Define spacecraft actors
+    sat1 = ActorBuilder.get_actor_scaffold("sat1", SpacecraftActor, pk.epoch(0))
+    sat2 = ActorBuilder.get_actor_scaffold("sat2", SpacecraftActor, pk.epoch(0))
+
+    # geostationary orbital parameters:
+    r = 6371000 + 35786000  # radius [km]
+    v = 3074.66  # velocity [m/s]
+
+    # To have a more symmetric case, let the actors be on same longitude as Earth magnetic dipole vector
+    longitude = -71.6 * np.pi / 180
+
+    # set orbits:
+    ActorBuilder.set_orbit(
+        sat1,
+        position=[
+            r * np.cos(np.pi / 2 + longitude),
+            r * np.sin(np.pi / 2 + longitude),
+            0,
+        ],
+        velocity=[
+            -v * np.sin(np.pi / 2 + longitude),
+            v * np.cos(np.pi / 2 + longitude),
+            0,
+        ],
+        epoch=pk.epoch(0),
+        central_body=earth,
+    )
+    ActorBuilder.set_orbit(
+        sat2,
+        position=[
+            r * np.cos(np.pi / 2 + longitude),
+            r * np.sin(np.pi / 2 + longitude),
+            0,
+        ],
+        velocity=[
+            -v * np.sin(np.pi / 2 + longitude),
+            v * np.cos(np.pi / 2 + longitude),
+            0,
+        ],
+        epoch=pk.epoch(0),
+        central_body=earth,
+    )
+    print([
+            r * np.cos(np.pi / 2 + longitude),
+            r * np.sin(np.pi / 2 + longitude),
+            0,
+        ],)
+    # set geometric model
+    ActorBuilder.set_geometric_model(sat1, mass=100)
+    ActorBuilder.set_geometric_model(sat2, mass=100)
+
+    # now, align the body magnetic dipole with the local Earth magnetic flux density vector
+    # Earth magnetic flux density vector at start position is approximately:
+
+    B = np.array([-3.18159529e-09, 1.02244882e-07, -3.72362170e-08])
+
+    B_direction = B / np.linalg.norm(B)
+    # define a very large dipole moment for magnetic actor to compensate for the low magnetic field at GEO orbit
+    m_body = 500  # Am²
+    actor_dipole = np.ndarray.tolist(B_direction * m_body)
+    initial_pointing_vector_body = np.ndarray.tolist(B_direction)
+
+    # set attitude models
+    ActorBuilder.set_attitude_model(
+        sat1,
+        actor_initial_angular_velocity=[0.0, 0.0, 0.0],
+        actor_pointing_vector_body=initial_pointing_vector_body,
+        actor_initial_attitude_in_rad=[0.0, 0.0, 0.0],
+        actor_residual_magnetic_field=actor_dipole,
+    )
+
+    ActorBuilder.set_attitude_model(
+        sat2,
+        actor_initial_angular_velocity=[0.0, 0.0, 0.0],
+        actor_pointing_vector_body=initial_pointing_vector_body,
+        actor_initial_attitude_in_rad=[0.0, 0.0, 0.0],
+        actor_residual_magnetic_field=[0.0, 0.0, 0.0],
+    )
+
+    # disturbances:
+    ActorBuilder.set_disturbances(sat1, magnetic=True)
+    ActorBuilder.set_disturbances(sat2, magnetic=True)
+
+    # Initial pointing vector in Earth inertial frame
+    initial_pointing_vector_eci = np.array(sat1.pointing_vector())
+    # Check if pointing vectors in Earth inertial frame are equal
+    assert np.all(sat1.pointing_vector() == sat2.pointing_vector())
+
+    # start simulation
+    sim = paseos.init_sim(sat1)
+    sim.add_known_actor(sat2)
+
+    for i in range(20):
+        sim.advance_time(200, 0)
+
+        # get Earth B vector at timestep
+        # Earth magnetic dipole moment:
+        m_earth = sat1._attitude_model.earth_magnetic_dipole_moment()
+        # parameters to calculate local B vector:
+        actor_position = np.array(sat1.get_position(sat1.local_time))
+        r = np.linalg.norm(actor_position)
+        r_hat = actor_position / r
+        # local B vector:
+        B = 1e-7 * (3 * np.dot(m_earth, r_hat) * r_hat - m_earth) / (r**3)
+
+        # B vector direction:
+        B_direction = B / np.linalg.norm(B)
+
+        # angle between the B vector and the actor's magnetic dipole vector (which is in the pointing vector direction):
+        angle = np.arccos(np.dot(B_direction, sat1.pointing_vector())) * 180/np.pi
+
+        # check if the magnetic actor dipole moment vector doesn't deviate more than 1° from the  B vector.
+        assert angle < 1
+
+    # check if the non-magnetic actor didn't rotate
+    assert np.all(sat2.pointing_vector() == initial_pointing_vector_eci)
+
+
+magnetic_disturbance_test()
diff --git a/test_attitude_code/test_attitude_plotting.py b/test_attitude_code/test_attitude_plotting.py
index 4168231..d6ef9b7 100644
--- a/test_attitude_code/test_attitude_plotting.py
+++ b/test_attitude_code/test_attitude_plotting.py
@@ -7,6 +7,7 @@
 import paseos
 from paseos.actors.spacecraft_actor import SpacecraftActor
 from paseos.actors.actor_builder import ActorBuilder
+from paseos.utils.reference_frame_transfer import rpy_to_body, eci_to_rpy
 
 matplotlib.use("Qt5Agg")
 
@@ -24,8 +25,9 @@
 R = 6371000 + 35786000
 v = 3074.66
 
-arr = np.array([-1.24217547e-09, 1.01735657e-07, -3.87201400e-08])
-initial_magn = np.ndarray.tolist(arr / np.linalg.norm(arr) * 20)
+#arr = np.array([-1.24217547e-09, 1.01735657e-07, -3.87201400e-08])
+arr = np.array([-3.18159529e-09, 1.02244882e-07, -3.72362170e-08])
+initial_magn = np.ndarray.tolist(arr / np.linalg.norm(arr) * 1000)
 initial_pointing = np.ndarray.tolist(arr / np.linalg.norm(arr))
 
 ActorBuilder.set_orbit(
@@ -46,7 +48,6 @@
 ActorBuilder.set_geometric_model(sat1, mass=100)
 ActorBuilder.set_geometric_model(sat2, mass=100)
 
-# when i = 21 in loop and advance time =100, pi/2000 rad/sec will rotate 180 deg about 1 axis
 ActorBuilder.set_attitude_model(
     sat1,
     actor_initial_angular_velocity=[0.0, 0.0, 0.0],
@@ -114,8 +115,18 @@
     # print("plotted attitude:", euler, " at position: ", pos, " pointing v: ", vector / 2e6)
 
     m = sat1._attitude_model.earth_magnetic_dipole_moment() * 6e-16
-    B = sat1._attitude_model._actor_magnetic_flux * 2e14
-    print(sat1._attitude_model._actor_magnetic_flux)
+
+    # get new Earth B vector
+    m_earth = sat1._attitude_model.earth_magnetic_dipole_moment()
+    actor_position = np.array(sat1.get_position(sat1.local_time))
+    actor_velocity = np.array(sat1.get_position_velocity(sat1.local_time)[1])
+    r = np.linalg.norm(actor_position)
+    r_hat = actor_position / r
+
+    B = 1e-7 * (3 * np.dot(m_earth, r_hat) * r_hat - m_earth) / (r**3)
+    #B = rpy_to_body(eci_to_rpy(B, actor_position, actor_velocity), sat1.attitude_in_rad())
+
+    B = B * 2e14
     # plot vectors
     # ax.quiver(pos[0], pos[1], pos[2], ang_vel[0], ang_vel[1], ang_vel[2], color="m")
     ax.quiver(pos[0], pos[1], pos[2], vector[0], vector[1], vector[2], linewidths=3, color="r")
@@ -125,8 +136,6 @@
     ax.quiver(pos[0], pos[1], pos[2], B[0], B[1], B[2], color="y")
 
     sim.advance_time(1000, 0)
-    # sim.advance_time(3446, 0)
-
 
 # 3D figure limits
 axmin = min([min(x), min(y), min(z)]) * 1.1

From ef25a18c0cddd8672c432606f5e70d460b8d41c6 Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Mon, 29 Jan 2024 17:59:33 +0100
Subject: [PATCH 09/16] test for magnetic disturbance finished

---
 paseos/tests/attitude_test.py | 95 +++++++++++++++++++++++------------
 1 file changed, 63 insertions(+), 32 deletions(-)

diff --git a/paseos/tests/attitude_test.py b/paseos/tests/attitude_test.py
index 37acb28..dd1982f 100644
--- a/paseos/tests/attitude_test.py
+++ b/paseos/tests/attitude_test.py
@@ -179,21 +179,57 @@ def magnetic_disturbance_test():
     magnetic field. One is non-magnetic and is expected to have a fixed attitude in the Earth inertial frame.
     The other (magnetic) actor should stay aligned with the Earth magnetic field.
     """
+
+    def flux_density_vector(actor, frame="eci"):
+        """Function to calculate the local magnetic field flux density vector (B) at the actor's location.
+        B vector is calculated multiple times. Use of this function for clarity.
+
+        Args:
+            actor (SpacecraftActor): Spacecraft actor
+            frame (string): B expressed in which frame (actor body frame or Earth-centered inertial frame)
+
+        Returns: B vector (np.ndarray)
+        """
+        # get Earth B vector at specific timestep
+
+        # Earth magnetic dipole moment:
+        m_earth = actor._attitude_model.earth_magnetic_dipole_moment()
+
+        # parameters to calculate local B vector:
+        actor_position = np.array(actor.get_position(actor.local_time))
+        r = np.linalg.norm(actor_position)
+        r_hat = actor_position / r
+
+        # local B vector:
+        B = 1e-7 * (3 * np.dot(m_earth, r_hat) * r_hat - m_earth) / (r**3)
+        if frame == "eci":
+            # return B vector in ECI frame
+            return B
+        elif frame == "body":
+            # transform B vector to the actor body frame and return
+            actor_velocity = np.array(actor.get_position_velocity(actor.local_time)[1])
+            actor_attitude = np.array(actor.attitude_in_rad())
+            return rpy_to_body(
+                eci_to_rpy(B, actor_position, actor_velocity), actor_attitude
+            )
+
     # Define central body
     earth = pk.planet.jpl_lp("earth")
 
     # Define spacecraft actors
+    # magnetic:
     sat1 = ActorBuilder.get_actor_scaffold("sat1", SpacecraftActor, pk.epoch(0))
+    # non-magnetic:
     sat2 = ActorBuilder.get_actor_scaffold("sat2", SpacecraftActor, pk.epoch(0))
 
-    # geostationary orbital parameters:
+    # Geostationary orbital parameters:
     r = 6371000 + 35786000  # radius [km]
     v = 3074.66  # velocity [m/s]
 
     # To have a more symmetric case, let the actors be on same longitude as Earth magnetic dipole vector
     longitude = -71.6 * np.pi / 180
 
-    # set orbits:
+    # Set orbits:
     ActorBuilder.set_orbit(
         sat1,
         position=[
@@ -224,33 +260,29 @@ def magnetic_disturbance_test():
         epoch=pk.epoch(0),
         central_body=earth,
     )
-    print([
-            r * np.cos(np.pi / 2 + longitude),
-            r * np.sin(np.pi / 2 + longitude),
-            0,
-        ],)
-    # set geometric model
+
+    # Set geometric model
     ActorBuilder.set_geometric_model(sat1, mass=100)
     ActorBuilder.set_geometric_model(sat2, mass=100)
 
-    # now, align the body magnetic dipole with the local Earth magnetic flux density vector
+    # 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])
 
-    B = np.array([-3.18159529e-09, 1.02244882e-07, -3.72362170e-08])
+    B_direction = B0 / np.linalg.norm(B0)
 
-    B_direction = B / np.linalg.norm(B)
-    # define a very large dipole moment for magnetic actor to compensate for the low magnetic field at GEO orbit
+    # Define a very large dipole moment for magnetic actor to compensate for the low magnetic field at GEO orbit
     m_body = 500  # Am²
     actor_dipole = np.ndarray.tolist(B_direction * m_body)
     initial_pointing_vector_body = np.ndarray.tolist(B_direction)
 
-    # set attitude models
+    # Set attitude models
     ActorBuilder.set_attitude_model(
         sat1,
         actor_initial_angular_velocity=[0.0, 0.0, 0.0],
         actor_pointing_vector_body=initial_pointing_vector_body,
         actor_initial_attitude_in_rad=[0.0, 0.0, 0.0],
-        actor_residual_magnetic_field=actor_dipole,
+        actor_residual_magnetic_field=actor_dipole,  # magnetic
     )
 
     ActorBuilder.set_attitude_model(
@@ -258,10 +290,10 @@ def magnetic_disturbance_test():
         actor_initial_angular_velocity=[0.0, 0.0, 0.0],
         actor_pointing_vector_body=initial_pointing_vector_body,
         actor_initial_attitude_in_rad=[0.0, 0.0, 0.0],
-        actor_residual_magnetic_field=[0.0, 0.0, 0.0],
+        actor_residual_magnetic_field=[0.0, 0.0, 0.0],  # non-magnetic
     )
 
-    # disturbances:
+    # Disturbances:
     ActorBuilder.set_disturbances(sat1, magnetic=True)
     ActorBuilder.set_disturbances(sat2, magnetic=True)
 
@@ -270,33 +302,32 @@ def magnetic_disturbance_test():
     # Check if pointing vectors in Earth inertial frame are equal
     assert np.all(sat1.pointing_vector() == sat2.pointing_vector())
 
-    # start simulation
+    # Start simulation
     sim = paseos.init_sim(sat1)
     sim.add_known_actor(sat2)
 
+    # Check if B0, used to define the actors' magnetic field orientations is the initial B vector orientation in sim.
+    assert np.all(np.isclose(B0, flux_density_vector(sat1, "body")))
+    assert np.all(np.isclose(B0, flux_density_vector(sat2, "body")))
+
+    # The magnetic actor will oscillate slightly, following the Earth's magnetic field lines. (check for multiple steps)
+    # The actor's magnetic field will not stay perfectly aligned with the Earth's field but needs to stay close.
     for i in range(20):
         sim.advance_time(200, 0)
-
-        # get Earth B vector at timestep
-        # Earth magnetic dipole moment:
-        m_earth = sat1._attitude_model.earth_magnetic_dipole_moment()
-        # parameters to calculate local B vector:
-        actor_position = np.array(sat1.get_position(sat1.local_time))
-        r = np.linalg.norm(actor_position)
-        r_hat = actor_position / r
-        # local B vector:
-        B = 1e-7 * (3 * np.dot(m_earth, r_hat) * r_hat - m_earth) / (r**3)
-
+        
+        # Get the magnetic flux density vector:
+        B = flux_density_vector(sat1)
+        
         # B vector direction:
         B_direction = B / np.linalg.norm(B)
 
-        # angle between the B vector and the actor's magnetic dipole vector (which is in the pointing vector direction):
-        angle = np.arccos(np.dot(B_direction, sat1.pointing_vector())) * 180/np.pi
+        # Angle between the B vector and the actor's magnetic dipole vector (which is in the pointing vector direction):
+        angle = np.arccos(np.dot(B_direction, sat1.pointing_vector())) * 180 / np.pi
 
-        # check if the magnetic actor dipole moment vector doesn't deviate more than 1° from the  B vector.
+        # Check if the magnetic actor dipole moment vector doesn't deviate more than 1° from the  B vector.
         assert angle < 1
 
-    # check if the non-magnetic actor didn't rotate
+    # Check if the non-magnetic actor didn't rotate
     assert np.all(sat2.pointing_vector() == initial_pointing_vector_eci)
 
 

From 7914692464ab4e414016f15928b338d38d1a9aee Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Mon, 29 Jan 2024 18:01:51 +0100
Subject: [PATCH 10/16] test for magnetic disturbance cleaned up

---
 paseos/tests/attitude_test.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/paseos/tests/attitude_test.py b/paseos/tests/attitude_test.py
index dd1982f..18b8425 100644
--- a/paseos/tests/attitude_test.py
+++ b/paseos/tests/attitude_test.py
@@ -174,15 +174,15 @@ def attitude_and_orbit_test():
 
 def magnetic_disturbance_test():
     """Tests the magnetic disturbance torques applied in the attitude model.
-    First: put two spacecraft actors in a geostationary orbit (disregarding the relative magnetic field rotation of the
-    Earth). Both actor's own magnetic dipole moment aligned with the local magnetic flux density vector of the Earth
+    Put two spacecraft actors in a geostationary orbit (disregarding the relative magnetic field rotation of the Earth).
+    Both actor's own magnetic dipole moment aligned with the local magnetic flux density vector of the Earth
     magnetic field. One is non-magnetic and is expected to have a fixed attitude in the Earth inertial frame.
     The other (magnetic) actor should stay aligned with the Earth magnetic field.
     """
 
     def flux_density_vector(actor, frame="eci"):
         """Function to calculate the local magnetic field flux density vector (B) at the actor's location.
-        B vector is calculated multiple times. Use of this function for clarity.
+        B vector is calculated multiple times. Use of this function for code clarity.
 
         Args:
             actor (SpacecraftActor): Spacecraft actor

From 21d963422f2a24e657b52338afc194ad4a013df9 Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Mon, 29 Jan 2024 18:26:25 +0100
Subject: [PATCH 11/16] deleted unnecessary lines

---
 paseos/tests/attitude_test.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/paseos/tests/attitude_test.py b/paseos/tests/attitude_test.py
index 18b8425..6c0cbff 100644
--- a/paseos/tests/attitude_test.py
+++ b/paseos/tests/attitude_test.py
@@ -151,7 +151,6 @@ def attitude_and_orbit_test():
         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])
@@ -329,6 +328,3 @@ def flux_density_vector(actor, frame="eci"):
 
     # Check if the non-magnetic actor didn't rotate
     assert np.all(sat2.pointing_vector() == initial_pointing_vector_eci)
-
-
-magnetic_disturbance_test()

From 8e37a3b8ba1e68d25646d144cef0322071f3ff55 Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Mon, 29 Jan 2024 19:15:08 +0100
Subject: [PATCH 12/16] no more next position in magnetic disturbance model.

---
 paseos/attitude/attitude_model.py | 14 +++-----------
 1 file changed, 3 insertions(+), 11 deletions(-)

diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index 2316c6a..a6c8d92 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -134,13 +134,6 @@ def calculate_disturbance_torque(self):
             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])
-        dt = 10
-        # time
-        t = self._actor.local_time
-        next_position = np.array(
-            self._actor.get_position(pk.epoch(t.mjd2000 + dt * pk.SEC2DAY, "mjd2000"))
-        )
-        position = np.array(self._actor.get_position(t))
 
         if self._actor.has_attitude_disturbances:
             if "aerodynamic" in self._actor.get_disturbances():
@@ -148,13 +141,12 @@ def calculate_disturbance_torque(self):
             if "gravitational" in self._actor.get_disturbances():
                 T += calculate_grav_torque()
             if "magnetic" in self._actor.get_disturbances():
+                time = self._actor.local_time
                 T += calculate_magnetic_torque(
                     m_earth=self.earth_magnetic_dipole_moment(),
                     m_sat=self._actor_residual_magnetic_field,
-                    position=next_position,
-                    velocity=self._actor.get_position_velocity(self._actor.local_time)[
-                        1
-                    ],
+                    position=self._actor.get_position(time),
+                    velocity=self._actor.get_position_velocity(time)[1],
                     attitude=self._actor_attitude_in_rad,
                 )
         return T

From f257760928e05b338687441c43dd271e09651960 Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Mon, 29 Jan 2024 23:25:43 +0100
Subject: [PATCH 13/16] commenting & cleaning up

---
 paseos/attitude/attitude_model.py            | 8 ++++++--
 test_attitude_code/test_attitude_plotting.py | 8 +++-----
 2 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index bc9d1eb..ca0fc09 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -59,6 +59,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
@@ -83,6 +85,8 @@ def __init__(
             np.array(self._actor.get_position(self._actor.local_time)),
             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):
@@ -94,7 +98,7 @@ def _nadir_vector(self):
         u = np.array(self._actor.get_position(self._actor.local_time))
         return -u / np.linalg.norm(u)
 
-    def earth_magnetic_dipole_moment(self):
+    def _earth_magnetic_dipole_moment(self):
         """Returns the Earth magnetic dipole moment vector 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.
@@ -142,7 +146,7 @@ def _calculate_disturbance_torque(self):
             if "magnetic" in self._actor.get_disturbances():
                 time = self._actor.local_time
                 T += calculate_magnetic_torque(
-                    m_earth=self.earth_magnetic_dipole_moment(),
+                    m_earth=self._earth_magnetic_dipole_moment(),
                     m_sat=self._actor_residual_magnetic_field,
                     position=self._actor.get_position(time),
                     velocity=self._actor.get_position_velocity(time)[1],
diff --git a/test_attitude_code/test_attitude_plotting.py b/test_attitude_code/test_attitude_plotting.py
index d6ef9b7..5d3bac1 100644
--- a/test_attitude_code/test_attitude_plotting.py
+++ b/test_attitude_code/test_attitude_plotting.py
@@ -7,7 +7,6 @@
 import paseos
 from paseos.actors.spacecraft_actor import SpacecraftActor
 from paseos.actors.actor_builder import ActorBuilder
-from paseos.utils.reference_frame_transfer import rpy_to_body, eci_to_rpy
 
 matplotlib.use("Qt5Agg")
 
@@ -104,7 +103,7 @@
     # angular velocity vector:
     # normalize first:
     ang_vel = sat1.angular_velocity()
-    if all(ang_vel == np.zeros(3)):
+    if np.all(ang_vel == np.zeros(3)):
         ang_vel = np.zeros(3)
     else:
         ang_vel = sat1.angular_velocity() / np.linalg.norm(sat1.angular_velocity())
@@ -114,17 +113,16 @@
 
     # print("plotted attitude:", euler, " at position: ", pos, " pointing v: ", vector / 2e6)
 
-    m = sat1._attitude_model.earth_magnetic_dipole_moment() * 6e-16
+    m = sat1._attitude_model._earth_magnetic_dipole_moment() * 6e-16
 
     # get new Earth B vector
-    m_earth = sat1._attitude_model.earth_magnetic_dipole_moment()
+    m_earth = sat1._attitude_model._earth_magnetic_dipole_moment()
     actor_position = np.array(sat1.get_position(sat1.local_time))
     actor_velocity = np.array(sat1.get_position_velocity(sat1.local_time)[1])
     r = np.linalg.norm(actor_position)
     r_hat = actor_position / r
 
     B = 1e-7 * (3 * np.dot(m_earth, r_hat) * r_hat - m_earth) / (r**3)
-    #B = rpy_to_body(eci_to_rpy(B, actor_position, actor_velocity), sat1.attitude_in_rad())
 
     B = B * 2e14
     # plot vectors

From e94520a7c7e92ec4eb332a3e991031251cfb382e Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Mon, 29 Jan 2024 23:45:42 +0100
Subject: [PATCH 14/16] todo added

---
 paseos/attitude/attitude_model.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index ca0fc09..3dd1423 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -139,6 +139,7 @@ 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():

From e51829e805d3218eb10a8bfbd19f52c45c9be515 Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Tue, 30 Jan 2024 00:33:33 +0100
Subject: [PATCH 15/16] Deleted file

---
 test_attitude_code/test_attitude_plotting.py | 152 -------------------
 1 file changed, 152 deletions(-)
 delete mode 100644 test_attitude_code/test_attitude_plotting.py

diff --git a/test_attitude_code/test_attitude_plotting.py b/test_attitude_code/test_attitude_plotting.py
deleted file mode 100644
index 5d3bac1..0000000
--- a/test_attitude_code/test_attitude_plotting.py
+++ /dev/null
@@ -1,152 +0,0 @@
-# We use pykep for orbit determination
-import pykep as pk
-import matplotlib
-import matplotlib.pyplot as plt
-import numpy as np
-
-import paseos
-from paseos.actors.spacecraft_actor import SpacecraftActor
-from paseos.actors.actor_builder import ActorBuilder
-
-matplotlib.use("Qt5Agg")
-
-# Define central body
-earth = pk.planet.jpl_lp("earth")
-
-# Define spacecraft actor
-sat1 = ActorBuilder.get_actor_scaffold("sat1", SpacecraftActor, pk.epoch(0))
-sat2 = ActorBuilder.get_actor_scaffold("sat2", SpacecraftActor, pk.epoch(0))
-# geostationary orbit:
-
-lat = 79.6 * np.pi / 180
-lon = -71.6 * np.pi / 180
-
-R = 6371000 + 35786000
-v = 3074.66
-
-#arr = np.array([-1.24217547e-09, 1.01735657e-07, -3.87201400e-08])
-arr = np.array([-3.18159529e-09, 1.02244882e-07, -3.72362170e-08])
-initial_magn = np.ndarray.tolist(arr / np.linalg.norm(arr) * 1000)
-initial_pointing = np.ndarray.tolist(arr / np.linalg.norm(arr))
-
-ActorBuilder.set_orbit(
-    sat1,
-    position=[R * np.cos(np.pi / 2 + lon), R * np.sin(np.pi / 2 + lon), 0],
-    velocity=[-v * np.sin(np.pi / 2 + lon), v * np.cos(np.pi / 2 + lon), 0],
-    epoch=pk.epoch(0),
-    central_body=earth,
-)
-ActorBuilder.set_orbit(
-    sat2,
-    position=[R * np.cos(np.pi / 2 + lon), R * np.sin(np.pi / 2 + lon), 0],
-    velocity=[-v * np.sin(np.pi / 2 + lon), v * np.cos(np.pi / 2 + lon), 0],
-    epoch=pk.epoch(0),
-    central_body=earth,
-)
-
-ActorBuilder.set_geometric_model(sat1, mass=100)
-ActorBuilder.set_geometric_model(sat2, mass=100)
-
-ActorBuilder.set_attitude_model(
-    sat1,
-    actor_initial_angular_velocity=[0.0, 0.0, 0.0],
-    actor_pointing_vector_body=initial_pointing,
-    actor_initial_attitude_in_rad=[0.0, 0.0, 0.0],
-    actor_residual_magnetic_field=initial_magn,
-)
-ActorBuilder.set_attitude_model(
-    sat2,
-    actor_initial_angular_velocity=[0.0, 0.0, 0.0],
-    actor_pointing_vector_body=initial_pointing,
-    actor_initial_attitude_in_rad=[0.0, 0.0, 0.0],
-    actor_residual_magnetic_field=[0.0, 0.0, 0.0],
-)
-
-# disturbances:
-ActorBuilder.set_disturbances(sat1, magnetic=True)
-ActorBuilder.set_disturbances(sat2, magnetic=True)
-
-
-sim = paseos.init_sim(sat1)
-sim.add_known_actor(sat2)
-plt.close()
-
-pos = []
-x = []
-y = []
-z = []
-pointing_vector = []
-
-fig = plt.figure()
-ax = fig.add_subplot(111, projection="3d")
-for i in range(46):
-    print("----------", i, "----------")
-    pos = sat1.get_position(sat1.local_time)
-    x.append(sat1.get_position(sat1.local_time)[0])
-    y.append(sat1.get_position(sat1.local_time)[1])
-    z.append(sat1.get_position(sat1.local_time)[2])
-
-    euler = sat1.attitude_in_rad()
-
-    # pointing vector:
-    vector = sat1.pointing_vector()
-    vector[np.isclose(vector, np.zeros(3))] = 0
-    # scale for plotting
-    vector = vector * 8e6
-
-    # pointing vector:
-    vector2 = sat2.pointing_vector()
-    vector2[np.isclose(vector2, np.zeros(3))] = 0
-    # scale for plotting
-    vector2 = vector2 * 8e6
-
-    # angular velocity vector:
-    # normalize first:
-    ang_vel = sat1.angular_velocity()
-    if np.all(ang_vel == np.zeros(3)):
-        ang_vel = np.zeros(3)
-    else:
-        ang_vel = sat1.angular_velocity() / np.linalg.norm(sat1.angular_velocity())
-        ang_vel[np.isclose(ang_vel, np.zeros(3))] = 0
-    # scale for plotting
-    ang_vel = ang_vel * 2e7
-
-    # print("plotted attitude:", euler, " at position: ", pos, " pointing v: ", vector / 2e6)
-
-    m = sat1._attitude_model._earth_magnetic_dipole_moment() * 6e-16
-
-    # get new Earth B vector
-    m_earth = sat1._attitude_model._earth_magnetic_dipole_moment()
-    actor_position = np.array(sat1.get_position(sat1.local_time))
-    actor_velocity = np.array(sat1.get_position_velocity(sat1.local_time)[1])
-    r = np.linalg.norm(actor_position)
-    r_hat = actor_position / r
-
-    B = 1e-7 * (3 * np.dot(m_earth, r_hat) * r_hat - m_earth) / (r**3)
-
-    B = B * 2e14
-    # plot vectors
-    # ax.quiver(pos[0], pos[1], pos[2], ang_vel[0], ang_vel[1], ang_vel[2], color="m")
-    ax.quiver(pos[0], pos[1], pos[2], vector[0], vector[1], vector[2], linewidths=3, color="r")
-    ax.quiver(pos[0], pos[1], pos[2], vector2[0], vector2[1], vector2[2], linewidths=3, color="m")
-    if not i % 10:
-        ax.quiver(0, 0, 0, m[0], m[1], m[2], color="y")
-    ax.quiver(pos[0], pos[1], pos[2], B[0], B[1], B[2], color="y")
-
-    sim.advance_time(1000, 0)
-
-# 3D figure limits
-axmin = min([min(x), min(y), min(z)]) * 1.1
-axmax = max([max(x), max(y), max(z)]) * 1.1
-
-ax.axes.set_xlim3d(left=axmin, right=axmax)
-ax.axes.set_ylim3d(bottom=axmin, top=axmax)
-ax.axes.set_zlim3d(bottom=axmin, top=axmax)
-
-ax.set_xlabel("x")
-ax.set_ylabel("y")
-ax.set_zlabel("z")
-
-ax.plot(x, y, z)
-ax.scatter(0, 0, 0)
-plt.show()

From 70a4e27daad98e15eb965c8e8df1e5da46250ea6 Mon Sep 17 00:00:00 2001
From: Mr-Medina <M.MedinaLeon@student.tudelft.nl>
Date: Tue, 30 Jan 2024 20:19:41 +0100
Subject: [PATCH 16/16] Moved magnetic_dipole_moment function to central body.

---
 paseos/attitude/attitude_model.py   | 35 +-----------
 paseos/central_body/central_body.py | 82 +++++++++++++++++++++++++++--
 paseos/tests/attitude_test.py       |  3 +-
 3 files changed, 79 insertions(+), 41 deletions(-)

diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index 3dd1423..170f6bf 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -1,6 +1,5 @@
 import numpy as np
 import pykep as pk
-from skyfield.api import wgs84, load
 
 
 from paseos.attitude.disturbance_calculations import (
@@ -98,38 +97,6 @@ def _nadir_vector(self):
         u = np.array(self._actor.get_position(self._actor.local_time))
         return -u / np.linalg.norm(u)
 
-    def _earth_magnetic_dipole_moment(self):
-        """Returns the Earth magnetic dipole moment vector 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.
-
-        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)
-
-        Returns: Time-dependent Earth dipole moment vector in ECI (np.ndarray): [mx, my, mz]
-        """
-        # North geomagnetic pole location
-        lat = 79.6
-        lon = -71.6
-
-        # Converting time to skyfield to use its API
-        t_skyfield = load.timescale().tt_jd(self._actor.local_time.jd)
-
-        # North geomagnetic pole location on Earth surface in cartesian coordinates
-        dipole_north_direction = np.array(
-            wgs84.latlon(lat, lon).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) * 7.79e22
-        )
-        return magnetic_dipole_moment
-
     def _calculate_disturbance_torque(self):
         """Compute total torque due to user specified disturbances.
 
@@ -147,7 +114,7 @@ def _calculate_disturbance_torque(self):
             if "magnetic" in self._actor.get_disturbances():
                 time = self._actor.local_time
                 T += calculate_magnetic_torque(
-                    m_earth=self._earth_magnetic_dipole_moment(),
+                    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],
diff --git a/paseos/central_body/central_body.py b/paseos/central_body/central_body.py
index d844525..60dab13 100644
--- a/paseos/central_body/central_body.py
+++ b/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
diff --git a/paseos/tests/attitude_test.py b/paseos/tests/attitude_test.py
index 2866d44..e7efcbf 100644
--- a/paseos/tests/attitude_test.py
+++ b/paseos/tests/attitude_test.py
@@ -186,9 +186,8 @@ def flux_density_vector(actor, frame="eci"):
         Returns: B vector (np.ndarray)
         """
         # get Earth B vector at specific timestep
-
         # Earth magnetic dipole moment:
-        m_earth = actor._attitude_model.earth_magnetic_dipole_moment()
+        m_earth = actor.central_body.magnetic_dipole_moment(actor.local_time)
 
         # parameters to calculate local B vector:
         actor_position = np.array(actor.get_position(actor.local_time))