From c792502d718ff562c842f9d7cc5f9253c9c9ab6b Mon Sep 17 00:00:00 2001
From: Moanalengkeek <69005558+Moanalengkeek@users.noreply.github.com>
Date: Mon, 22 Jan 2024 12:55:14 +0100
Subject: [PATCH 1/6] Created disturbance model function and partly implemented
 it in the attitude model

---
 paseos/attitude/attitude_model.py           | 16 ++++++++--
 paseos/attitude/disturbance_calculations.py | 35 ++++++++++++++++++---
 2 files changed, 44 insertions(+), 7 deletions(-)

diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index 17ac692..b9211b3 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -12,6 +12,7 @@
     rodrigues_rotation,
     get_rpy_angles,
     rotate_body_vectors,
+    rpy_to_body,
 )
 
 
@@ -80,18 +81,29 @@ def nadir_vector(self):
         u = np.array(self._actor.get_position(self._actor.local_time))
         return -u / np.linalg.norm(u)
 
-    def calculate_disturbance_torque(self):
+    def calculate_disturbance_torque(self, position, velocity, euler_angles):
         """Compute total torque due to user specified disturbances
 
+        Args:
+            position (np.ndarray): position vector of RPY reference frame wrt ECI frame
+            velocity (np.ndarray): velocity of the spacecraft in earth reference frame, centered on spacecraft
+            euler_angles (np.ndarray): [roll, pitch, yaw] in radians
+
         Returns:
             np.array([Tx, Ty, Tz]): total combined torques in Nm expressed in the spacecraft body frame
         """
 
+        # Transform the earth rotation vector to the body reference frame, assuming the rotation vector is the z-axis
+        # of the earth-centered-inertial (eci) frame
+        earth_rotation_vector_in_rpy = eci_to_rpy(np.array([0, 0, 1]), position, velocity)
+        earth_rotation_vector_in_body = rpy_to_body(earth_rotation_vector_in_rpy, euler_angles)
+
         T = np.array([0.0, 0.0, 0.0])
         if "aerodynamic" in self._actor.get_disturbances():
             T += calculate_aero_torque()
         if "gravitational" in self._actor.get_disturbances():
-            T += calculate_grav_torque()
+            T += calculate_grav_torque(self._actor_pointing_vector_body,earth_rotation_vector_in_body,
+                                       self._actor._moment_of_inertia, self._actor._previous_altitude)
         if "magnetic" in self._actor.get_disturbances():
             T += calculate_magnetic_torque()
         return T
diff --git a/paseos/attitude/disturbance_calculations.py b/paseos/attitude/disturbance_calculations.py
index 43f6fe2..cb3d092 100644
--- a/paseos/attitude/disturbance_calculations.py
+++ b/paseos/attitude/disturbance_calculations.py
@@ -12,11 +12,36 @@ def calculate_aero_torque():
     return np.array(T)
 
 
-def calculate_grav_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_grav_torque(u_r, u_n, J, h):
+    """
+    Equation for gravity gradient torque with up to J2 effect from:
+    https://doi.org/10.1016/j.asr.2018.06.025, chapter 3.3
+    This function currently only works for Earth centered orbits
+
+    Args:
+        u_r (np.array): Unit vector pointing from Satellite center of gravity to Earth's center of gravity
+        u_n (np.array): Unit vector along the Earth's rotation axis, in the spacecraft body frame
+        J (np.array): The satellites moment of inertia, in the form of [[Ixx Ixy Ixz]
+                                                                        [Iyx Iyy Iyx]
+                                                                        [Izx Izy Izz]]
+        h (float): The altitude of the spacecraft in m
+
+    Returns:
+        np.array: total gravitational torques in Nm expressed in the spacecraft body frame
+    """
+    # Constants
+    mu = 3.986004418e14  # Earth's gravitational parameter, [m^3/s^2]
+    J2 = 1.0826267e-3  # Earth's J2 coefficient, from https://ocw.tudelft.nl/wp-content/uploads/AE2104-Orbital-Mechanics-Slides_8.pdf
+    Re = 6371000  # Earth's radius, [m]
+
+    # Simulation parameters
+    r = h + Re  # Radial distance from Satellite center of gravity to Earth's center of gravity, [m]
+
+    tg_term_1 = np.cross((3 * mu * u_r / (r ** 3)), J * u_r)
+    tg_term_2 = 30 * np.dot(u_r, u_n) * (np.cross(u_n, J * u_r) + np.cross(u_r, J * u_n))
+    tg_term_3 = np.cross((15 - 105 * np.dot(u_r, u_n) ** 2 * u_r), J * u_r) + np.cross(6 * u_n, J * u_r)
+    tg = tg_term_1 + mu * J2 * Re ** 2 / (2 * r ** 5) * (tg_term_2 + tg_term_3)
+    return np.array(tg)
 
 
 def calculate_magnetic_torque():

From acb6aeeaee0d372dc02f64ad8bb8ab33df48625f Mon Sep 17 00:00:00 2001
From: Moanalengkeek <69005558+Moanalengkeek@users.noreply.github.com>
Date: Thu, 25 Jan 2024 15:57:15 +0100
Subject: [PATCH 2/6] Updates gravity model implementation

---
 paseos/attitude/attitude_model.py           | 6 ++++--
 paseos/attitude/disturbance_calculations.py | 3 ++-
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index b9211b3..4c4610f 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -103,7 +103,7 @@ def calculate_disturbance_torque(self, position, velocity, euler_angles):
             T += calculate_aero_torque()
         if "gravitational" in self._actor.get_disturbances():
             T += calculate_grav_torque(self._actor_pointing_vector_body,earth_rotation_vector_in_body,
-                                       self._actor._moment_of_inertia, self._actor._previous_altitude)
+                                       self._actor._moment_of_inertia, position[0])
         if "magnetic" in self._actor.get_disturbances():
             T += calculate_magnetic_torque()
         return T
@@ -118,7 +118,9 @@ def calculate_angular_acceleration(self):
         # 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
         self._actor_angular_acceleration = np.linalg.inv(I) @ (
-            self.calculate_disturbance_torque()
+            self.calculate_disturbance_torque(position=np.array(self._actor.get_position(self._actor.local_time)),
+                                              velocity=np.array(self._actor.get_position_velocity(self._actor.local_time)[1]),
+                                              euler_angles=self._actor_attitude_in_rad)
             - np.cross(self._actor_angular_velocity, I @ self._actor_angular_velocity)
         )
 
diff --git a/paseos/attitude/disturbance_calculations.py b/paseos/attitude/disturbance_calculations.py
index cb3d092..4bdc9fe 100644
--- a/paseos/attitude/disturbance_calculations.py
+++ b/paseos/attitude/disturbance_calculations.py
@@ -41,7 +41,8 @@ def calculate_grav_torque(u_r, u_n, J, h):
     tg_term_2 = 30 * np.dot(u_r, u_n) * (np.cross(u_n, J * u_r) + np.cross(u_r, J * u_n))
     tg_term_3 = np.cross((15 - 105 * np.dot(u_r, u_n) ** 2 * u_r), J * u_r) + np.cross(6 * u_n, J * u_r)
     tg = tg_term_1 + mu * J2 * Re ** 2 / (2 * r ** 5) * (tg_term_2 + tg_term_3)
-    return np.array(tg)
+    T = [tg[0,0], tg[1,1], tg[2,2]]
+    return np.array(T)
 
 
 def calculate_magnetic_torque():

From d31310fe541bce9caee0b7dd17666349d76f950f Mon Sep 17 00:00:00 2001
From: Moanalengkeek <69005558+Moanalengkeek@users.noreply.github.com>
Date: Sat, 27 Jan 2024 16:08:40 +0100
Subject: [PATCH 3/6] Added tests and removed bugs from gravity disturbance
 model

---
 paseos/actors/actor_builder.py              |  1 +
 paseos/attitude/attitude_model.py           |  6 +-
 paseos/attitude/disturbance_calculations.py | 15 ++--
 paseos/tests/attitude_test.py               | 92 +++++++++++++++++++++
 4 files changed, 103 insertions(+), 11 deletions(-)
 create mode 100644 paseos/tests/attitude_test.py

diff --git a/paseos/actors/actor_builder.py b/paseos/actors/actor_builder.py
index 350ca88..51a6424 100644
--- a/paseos/actors/actor_builder.py
+++ b/paseos/actors/actor_builder.py
@@ -333,6 +333,7 @@ def set_geometric_model(
             local_actor=actor, actor_mass=mass, vertices=vertices, faces=faces, scale=scale
         )
         actor._mesh = geometric_model.set_mesh()
+        actor._center_of_gravity = geometric_model.find_center_of_gravity()
         actor._moment_of_inertia = geometric_model.find_moment_of_inertia
 
     @staticmethod
diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index 4c4610f..91d0dad 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -102,8 +102,10 @@ def calculate_disturbance_torque(self, position, velocity, euler_angles):
         if "aerodynamic" in self._actor.get_disturbances():
             T += calculate_aero_torque()
         if "gravitational" in self._actor.get_disturbances():
-            T += calculate_grav_torque(self._actor_pointing_vector_body,earth_rotation_vector_in_body,
-                                       self._actor._moment_of_inertia, position[0])
+            nadir_vector_in_rpy = eci_to_rpy(self.nadir_vector(), position, velocity)
+            nadir_vector_in_body = rpy_to_body(nadir_vector_in_rpy, euler_angles)
+            T += calculate_grav_torque(nadir_vector_in_body,earth_rotation_vector_in_body,
+                                       self._actor._moment_of_inertia, np.linalg.norm(position))
         if "magnetic" in self._actor.get_disturbances():
             T += calculate_magnetic_torque()
         return T
diff --git a/paseos/attitude/disturbance_calculations.py b/paseos/attitude/disturbance_calculations.py
index 4bdc9fe..3c1b6e2 100644
--- a/paseos/attitude/disturbance_calculations.py
+++ b/paseos/attitude/disturbance_calculations.py
@@ -12,7 +12,7 @@ def calculate_aero_torque():
     return np.array(T)
 
 
-def calculate_grav_torque(u_r, u_n, J, h):
+def calculate_grav_torque(u_r, u_n, J, r):
     """
     Equation for gravity gradient torque with up to J2 effect from:
     https://doi.org/10.1016/j.asr.2018.06.025, chapter 3.3
@@ -24,7 +24,7 @@ def calculate_grav_torque(u_r, u_n, J, h):
         J (np.array): The satellites moment of inertia, in the form of [[Ixx Ixy Ixz]
                                                                         [Iyx Iyy Iyx]
                                                                         [Izx Izy Izz]]
-        h (float): The altitude of the spacecraft in m
+        h (float): The distance from the center of the Earth to the satellite
 
     Returns:
         np.array: total gravitational torques in Nm expressed in the spacecraft body frame
@@ -34,15 +34,12 @@ def calculate_grav_torque(u_r, u_n, J, h):
     J2 = 1.0826267e-3  # Earth's J2 coefficient, from https://ocw.tudelft.nl/wp-content/uploads/AE2104-Orbital-Mechanics-Slides_8.pdf
     Re = 6371000  # Earth's radius, [m]
 
-    # Simulation parameters
-    r = h + Re  # Radial distance from Satellite center of gravity to Earth's center of gravity, [m]
 
-    tg_term_1 = np.cross((3 * mu * u_r / (r ** 3)), J * u_r)
-    tg_term_2 = 30 * np.dot(u_r, u_n) * (np.cross(u_n, J * u_r) + np.cross(u_r, J * u_n))
-    tg_term_3 = np.cross((15 - 105 * np.dot(u_r, u_n) ** 2 * u_r), J * u_r) + np.cross(6 * u_n, J * u_r)
+    tg_term_1 = (3 * mu / (r ** 3))*np.cross(u_r, np.dot(J,u_r))
+    tg_term_2 = 30 * np.dot(u_r, u_n)*(np.cross(u_n, np.dot(J,u_r)) + np.cross(u_r, np.dot(J,u_n)))
+    tg_term_3 = np.cross((15 - 105 * np.dot(u_r, u_n) ** 2) * u_r, np.dot(J,u_r)) + np.cross(6 * u_n, np.dot(J ,u_n))
     tg = tg_term_1 + mu * J2 * Re ** 2 / (2 * r ** 5) * (tg_term_2 + tg_term_3)
-    T = [tg[0,0], tg[1,1], tg[2,2]]
-    return np.array(T)
+    return np.array(tg)
 
 
 def calculate_magnetic_torque():
diff --git a/paseos/tests/attitude_test.py b/paseos/tests/attitude_test.py
new file mode 100644
index 0000000..8a4b3d0
--- /dev/null
+++ b/paseos/tests/attitude_test.py
@@ -0,0 +1,92 @@
+import numpy as np
+import pykep as pk
+import sys
+
+sys.path.append("../..")
+import paseos
+from paseos import ActorBuilder, SpacecraftActor, load_default_cfg
+
+
+def gravity_disturbance_cube_test():
+    """This test checks whether a 3-axis symmetric, uniform body (a cube with constant density, and cg at origin)
+    creates no angular acceleration/velocity due to gravity"""
+    earth = pk.planet.jpl_lp("earth")
+
+    # Define local actor
+    sat1 = ActorBuilder.get_actor_scaffold("sat1", SpacecraftActor, pk.epoch(0))
+    ActorBuilder.set_orbit(sat1, [7000000, 0, 0], [0, 8000.0, 0], pk.epoch(0), earth)
+    ActorBuilder.set_geometric_model(sat1, mass=100)
+    ActorBuilder.set_attitude_model(sat1)
+    ActorBuilder.set_disturbances(sat1, gravitational=True)
+
+    # init simulation
+    sim = paseos.init_sim(sat1)
+
+    # Check initial conditions
+    assert np.all(sat1._attitude_model._actor_angular_velocity == 0.0)
+
+    # run simulation for 1 period
+    orbital_period = 2 * np.pi * np.sqrt((6371000 + 7000000) ** 3 / 3.986004418e14)
+    sim.advance_time(orbital_period, 0)
+    nadir = sat1._attitude_model.nadir_vector()
+
+    # check conditions after 1 orbit
+    assert np.all(sat1._attitude_model._actor_angular_acceleration == 0.0)
+
+def gravity_disturbance_pole_test():
+    """This test checks whether a 2-axis symmetric, uniform body (a pole (10x1x1) with constant density, and cg at
+    origin) stabilises in orbit due to gravitational acceleration
+    It additionally checks the implementation of custom meshes of the geometric model"""
+
+    vertices = [
+                [-5, -0.5, -0.5],
+                [-5, -0.5, 0.5],
+        [-5, 0.5, -0.5],
+        [-5, 0.5, 0.5],
+        [5, -0.5, -0.5],
+        [5, -0.5, 0.5],
+        [5, 0.5, -0.5],
+        [5, 0.5, 0.5],
+    ]
+    faces = [
+        [0, 1, 3],
+        [0, 3, 2],
+        [0, 2, 6],
+        [0, 6, 4],
+        [1, 5, 3],
+        [3, 5, 7],
+        [2, 3, 7],
+        [2, 7, 6],
+        [4, 6, 7],
+        [4, 7, 5],
+        [0, 4, 1],
+        [1, 4, 5],
+    ]
+
+    earth = pk.planet.jpl_lp("earth")
+
+    # Define local actor
+    sat1 = ActorBuilder.get_actor_scaffold("sat1", SpacecraftActor, pk.epoch(0))
+    ActorBuilder.set_orbit(sat1, [7000000, 0, 0], [0, 8000.0, 0], pk.epoch(0), earth)
+    ActorBuilder.set_geometric_model(sat1, mass=100, vertices=vertices, faces=faces)
+    orbital_period = 2 * np.pi * np.sqrt((6371000 + 7000000) ** 3 / 3.986004418e14)
+    ActorBuilder.set_attitude_model(sat1)#, actor_initial_angular_velocity=[0,2*np.pi/orbital_period,0])
+    ActorBuilder.set_disturbances(sat1, gravitational=True)
+
+    # init simulation
+    cfg = load_default_cfg()  # loading cfg to modify defaults
+    cfg.sim.dt = 100.0  # setting higher timestep to run things quickly
+    sim = paseos.init_sim(sat1, cfg)
+
+
+    # Check initial conditions
+    assert np.all(sat1._attitude_model._actor_attitude_in_rad == 0.0)
+
+    # run simulation for 1 period
+    for i in range(11):
+        sim.advance_time(orbital_period*0.1, 0)
+
+    # check conditions after 0.1 orbit, satellite should have acceleration around y-axis to align pole towards earth
+    assert np.round(sat1._attitude_model._actor_angular_acceleration[0],10) == 0.0
+    assert not np.round(sat1._attitude_model._actor_angular_acceleration[1],10) == 0.0
+

From a2cb14164602b22b1b6d8cd2ae55129447ef5da3 Mon Sep 17 00:00:00 2001
From: Moanalengkeek <69005558+Moanalengkeek@users.noreply.github.com>
Date: Tue, 30 Jan 2024 10:46:35 +0100
Subject: [PATCH 4/6] Merged final attitude model inot this branch, and
 implemented changes in testing to account for that

---
 paseos/attitude/attitude_model.py | 7 ++++---
 paseos/tests/attitude_test.py     | 6 ++++--
 2 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index 9de7749..efba0c7 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -106,8 +106,7 @@ def calculate_disturbance_torque(self, position, velocity, euler_angles):
 
         # Transform the earth rotation vector to the body reference frame, assuming the rotation vector is the z-axis
         # of the earth-centered-inertial (eci) frame
-        earth_rotation_vector_in_rpy = eci_to_rpy(np.array([0, 0, 1]), position, velocity)
-        earth_rotation_vector_in_body = rpy_to_body(earth_rotation_vector_in_rpy, euler_angles)
+
 
         T = np.array([0.0, 0.0, 0.0])
 
@@ -117,8 +116,10 @@ def calculate_disturbance_torque(self, position, velocity, euler_angles):
             if "gravitational" in self._actor.get_disturbances():
                 nadir_vector_in_rpy = eci_to_rpy(self.nadir_vector(), position, velocity)
                 nadir_vector_in_body = rpy_to_body(nadir_vector_in_rpy, euler_angles)
+                earth_rotation_vector_in_rpy = eci_to_rpy(np.array([0, 0, 1]), position, velocity)
+                earth_rotation_vector_in_body = rpy_to_body(earth_rotation_vector_in_rpy, euler_angles)
                 T += calculate_grav_torque(nadir_vector_in_body,earth_rotation_vector_in_body,
-                                           self._actor._moment_of_inertia, np.linalg.norm(position))
+                                           self._actor._moment_of_inertia(), np.linalg.norm(position))
             if "magnetic" in self._actor.get_disturbances():
                 T += calculate_magnetic_torque()
         return T
diff --git a/paseos/tests/attitude_test.py b/paseos/tests/attitude_test.py
index a0dafc7..33b64b1 100644
--- a/paseos/tests/attitude_test.py
+++ b/paseos/tests/attitude_test.py
@@ -32,7 +32,8 @@ def gravity_disturbance_cube_test():
     nadir = sat1._attitude_model.nadir_vector()
 
     # check conditions after 1 orbit
-    assert np.all(sat1._attitude_model._actor_angular_acceleration == 0.0)
+    assert np.all(np.round(sat1._attitude_model._actor_angular_acceleration,10) == 0.0)
+
 
 def gravity_disturbance_pole_test():
     """This test checks whether a 2-axis symmetric, uniform body (a pole (10x1x1) with constant density, and cg at
@@ -91,6 +92,7 @@ def gravity_disturbance_pole_test():
     assert np.round(sat1._attitude_model._actor_angular_acceleration[0],10) == 0.0
     assert not np.round(sat1._attitude_model._actor_angular_acceleration[1],10) == 0.0
 
+
 def attitude_model_test():
     """Testing the attitude model with no disturbances and known angular velocity.
     One actor has orbit in Earth inertial x-y plane (equatorial) with initial velocity which rotates the actor with 180°
@@ -246,4 +248,4 @@ def attitude_and_orbit_test():
 
     # Testing the simulation went as intended
     assert sat1._attitude_model._actor_pointing_vector_body[2] == 1.0
-    assert vector[0] == -1.0
\ No newline at end of file
+    assert vector[0] == -1.0

From 3f1319071e09dda26cee9c9e2c1d4cdc2a22669f Mon Sep 17 00:00:00 2001
From: Moanalengkeek <69005558+Moanalengkeek@users.noreply.github.com>
Date: Wed, 31 Jan 2024 16:51:16 +0100
Subject: [PATCH 5/6] Update paseos/attitude/attitude_model.py

Co-authored-by: Gabriele Meoni <70584239+GabrieleMeoni@users.noreply.github.com>
---
 paseos/attitude/attitude_model.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index 82265c3..6d8367b 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -93,7 +93,7 @@ def _nadir_vector(self):
 
 
     def calculate_disturbance_torque(self, position, velocity, euler_angles):
-        """Compute total torque due to user specified disturbances
+        """Compute total torque due to user-specified disturbances.
 
         Args:
             position (np.ndarray): position vector of RPY reference frame wrt ECI frame

From a9dc4ab46697a16371ebd10b0ad4d383644342ec Mon Sep 17 00:00:00 2001
From: Moanalengkeek <69005558+Moanalengkeek@users.noreply.github.com>
Date: Wed, 31 Jan 2024 17:01:44 +0100
Subject: [PATCH 6/6] Final commits

---
 paseos/attitude/attitude_model.py           | 2 +-
 paseos/attitude/disturbance_calculations.py | 8 ++++----
 paseos/tests/attitude_test.py               | 6 ++++--
 3 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/paseos/attitude/attitude_model.py b/paseos/attitude/attitude_model.py
index efba0c7..3a63332 100644
--- a/paseos/attitude/attitude_model.py
+++ b/paseos/attitude/attitude_model.py
@@ -118,7 +118,7 @@ def calculate_disturbance_torque(self, position, velocity, euler_angles):
                 nadir_vector_in_body = rpy_to_body(nadir_vector_in_rpy, euler_angles)
                 earth_rotation_vector_in_rpy = eci_to_rpy(np.array([0, 0, 1]), position, velocity)
                 earth_rotation_vector_in_body = rpy_to_body(earth_rotation_vector_in_rpy, euler_angles)
-                T += calculate_grav_torque(nadir_vector_in_body,earth_rotation_vector_in_body,
+                T += calculate_grav_torque(self._actor.central_body.planet, nadir_vector_in_body,earth_rotation_vector_in_body,
                                            self._actor._moment_of_inertia(), np.linalg.norm(position))
             if "magnetic" in self._actor.get_disturbances():
                 T += calculate_magnetic_torque()
diff --git a/paseos/attitude/disturbance_calculations.py b/paseos/attitude/disturbance_calculations.py
index 3c1b6e2..8228271 100644
--- a/paseos/attitude/disturbance_calculations.py
+++ b/paseos/attitude/disturbance_calculations.py
@@ -12,7 +12,7 @@ def calculate_aero_torque():
     return np.array(T)
 
 
-def calculate_grav_torque(u_r, u_n, J, r):
+def calculate_grav_torque(central_body, u_r, u_n, J, r):
     """
     Equation for gravity gradient torque with up to J2 effect from:
     https://doi.org/10.1016/j.asr.2018.06.025, chapter 3.3
@@ -24,15 +24,15 @@ def calculate_grav_torque(u_r, u_n, J, r):
         J (np.array): The satellites moment of inertia, in the form of [[Ixx Ixy Ixz]
                                                                         [Iyx Iyy Iyx]
                                                                         [Izx Izy Izz]]
-        h (float): The distance from the center of the Earth to the satellite
+        r (float): The distance from the center of the Earth to the satellite
 
     Returns:
         np.array: total gravitational torques in Nm expressed in the spacecraft body frame
     """
     # Constants
-    mu = 3.986004418e14  # Earth's gravitational parameter, [m^3/s^2]
+    mu = central_body.mu_self  # Earth's gravitational parameter, [m^3/s^2]
     J2 = 1.0826267e-3  # Earth's J2 coefficient, from https://ocw.tudelft.nl/wp-content/uploads/AE2104-Orbital-Mechanics-Slides_8.pdf
-    Re = 6371000  # Earth's radius, [m]
+    Re = central_body.radius  # Earth's radius, [m]
 
 
     tg_term_1 = (3 * mu / (r ** 3))*np.cross(u_r, np.dot(J,u_r))
diff --git a/paseos/tests/attitude_test.py b/paseos/tests/attitude_test.py
index 33b64b1..18d136d 100644
--- a/paseos/tests/attitude_test.py
+++ b/paseos/tests/attitude_test.py
@@ -10,7 +10,8 @@
 
 def gravity_disturbance_cube_test():
     """This test checks whether a 3-axis symmetric, uniform body (a cube with constant density, and cg at origin)
-    creates no angular acceleration/velocity due to gravity"""
+    creates no angular acceleration/velocity due to gravity. The spacecraft is initially positioned with the z-axis
+    aligned with the nadir vector."""
     earth = pk.planet.jpl_lp("earth")
 
     # Define local actor
@@ -37,7 +38,8 @@ def gravity_disturbance_cube_test():
 
 def gravity_disturbance_pole_test():
     """This test checks whether a 2-axis symmetric, uniform body (a pole (10x1x1) with constant density, and cg at
-    origin) stabilises in orbit due to gravitational acceleration
+    origin) stabilises in orbit due to gravitational acceleration. The attitude at time zero is the z-axis pointing
+    towards earth. Hence, the accelerations should occur only in the y-axis.
     It additionally checks the implementation of custom meshes of the geometric model"""
 
     vertices = [