Moved magnetic_dipole_moment function to central body.

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],

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

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))