diff --git a/jobs/src/integrators.py b/jobs/src/integrators.py
index 79f0ec7b25567da4f0f182b81792e5a835564102..d0d786844dd8e2e534da0dbd330665ec9c7fa6a5 100644
--- a/jobs/src/integrators.py
+++ b/jobs/src/integrators.py
@@ -1,5 +1,5 @@
 """
-This module contains 2 integrators, which is at least 
+This module contains 1 integrator, which is at least 
 of second consistency order and approrimately preserves 
 the energy over long times.
 """
diff --git a/jobs/tests/test_bht_algorithm_2D.py b/jobs/tests/test_bht_algorithm_2D.py
index 3ea2ec12b424b1bb60932e0695b0b2b7af966699..3a409f6e61627b3b11bfae72267abed230ca4a09 100644
--- a/jobs/tests/test_bht_algorithm_2D.py
+++ b/jobs/tests/test_bht_algorithm_2D.py
@@ -74,7 +74,7 @@ class IntegratorTest(unittest.TestCase):
             0,
             R_SE * 2 * np.pi / T_E,
             0,
-            1 / M_E * M_E * R_SE * 2 * np.pi / T_E + 1 * R_L * 2 * np.pi / T_L,
+            R_SE * 2 * np.pi / T_E + 1 * R_L * 2 * np.pi / T_L,
         ])
 
         system = GravitationalSystem(r0=x0[:4],
@@ -89,11 +89,15 @@ class IntegratorTest(unittest.TestCase):
         ## plotting trajectories
         plt.figure()
         plt.plot(q[:, 0], q[:, 1], label="Earth")
+        plt.xlabel("x/AU")
+        plt.ylabel("y/AU")
         plt.legend()
         plt.show()
 
         plt.figure()
         plt.plot(q[:, 2] - q[:, 0], q[:, 3] - q[:, 1], label="Moon seen from Earth")
+        plt.xlabel("x/AU")
+        plt.ylabel("y/AU")
         plt.legend()
         plt.show()
 
@@ -113,8 +117,9 @@ class IntegratorTest(unittest.TestCase):
         self.assertTrue(np.greater(1e-8 + np.zeros(L.shape[0]), L - L[0]).all())
 
         ## checking error
-        dts = [dt, 2 * dt, 4 * dt]
-        errors = []
+        dts = [dt, 2 * dt, 4 * dt, 10 * dt, 50 * dt]
+        errors_E = []
+        errors_P = []
         ts = []
         for i in dts:
             system = GravitationalSystem(r0=x0[:4],
@@ -131,16 +136,36 @@ class IntegratorTest(unittest.TestCase):
             -G * M_S * M_E / np.linalg.norm(q_t[:,:2], axis=1) - G * M_S * M_L / np.linalg.norm(q_t[:,2:], axis=1) + \
             -G * M_E * M_L / np.linalg.norm(q_t[:,2:] - q_t[:,:2], axis=1)
 
-            errors.append((H - H[0]) / i**2)
+            P = p_t[:, :2] + p_t[:, 2:]
+
+            errors_E.append((H - H[0]) / i**2)
+            errors_P.append(np.linalg.norm(P, axis=1) - np.linalg.norm(P[0]))
             ts.append(t)
 
         plt.figure()
-        plt.plot(ts[0], errors[0], label="dt")
-        plt.plot(ts[1], errors[1], linestyle='--', label="2*dt")
-        plt.plot(ts[2], errors[2], linestyle=':', label="4*dt")
+        plt.plot(ts[0], errors_E[0], label="dt")
+        plt.plot(ts[1], errors_E[1], linestyle='--', label="2*dt")
+        plt.plot(ts[2], errors_E[2], linestyle=':', label="4*dt")
+        plt.plot(ts[3], errors_E[3], linestyle='-.', label="10*dt")
+        plt.plot(ts[4], errors_E[4], linestyle=':', label="50*dt")
         plt.xlabel("$t$")
         plt.ylabel("$\delta E(t)/(\Delta t)^2$")
         plt.legend()
+        plt.title("Energy error")
+        plt.tight_layout()
+        plt.show()
+
+        plt.figure()
+        plt.plot(ts[0], errors_P[0], label="dt")
+        plt.plot(ts[1], errors_P[1], linestyle='--', label="2*dt")
+        plt.plot(ts[2], errors_P[2], linestyle=':', label="4*dt")
+        plt.plot(ts[3], errors_P[3], linestyle='-.', label="10*dt")
+        plt.plot(ts[4], errors_P[4], linestyle=':', label="50*dt")
+        plt.xlabel("$t$")
+        plt.ylabel("$P(t)-P(0)$")
+        plt.legend()
+        plt.title("Linear momentum error")
+        plt.tight_layout()
         plt.show()
 
         ## checking time reversal: p -> -p
@@ -156,6 +181,29 @@ class IntegratorTest(unittest.TestCase):
         self.assertTrue(np.greater(1e-10 + np.zeros(4), q_reverse[-1] - q[0]).all())
         self.assertTrue(np.greater(1e-10 + np.zeros(4), p_reverse[-1] - p[0]).all())
 
+        ## plotting trajectories
+        plt.figure()
+        plt.plot(q[:, 0], q[:, 1], label="Forward")
+        plt.plot(q_reverse[:, 0], q_reverse[:, 1], '--', label="Backward")
+
+        plt.xlabel("x/AU")
+        plt.ylabel("y/AU")
+        plt.title("Earth")
+        plt.tight_layout()
+        plt.legend()
+        plt.show()
+
+        plt.figure()
+        plt.plot(q[:, 2] - q[:, 0], q[:, 3] - q[:, 1], label="Forward")
+        plt.plot(q_reverse[:, 2] - q_reverse[:, 0], q_reverse[:, 3] - q_reverse[:, 1], '--', label="Backward")
+
+        plt.xlabel("x/AU")
+        plt.ylabel("y/AU")
+        plt.title("Moon")
+        plt.tight_layout()
+        plt.legend()
+        plt.show()
+
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/jobs/tests/test_integrator.py b/jobs/tests/test_integrator.py
index e9bb8894eeaad42d7bd933e9396c836806f67d86..383182e3155e8b79af26703903e52434674ee6ee 100644
--- a/jobs/tests/test_integrator.py
+++ b/jobs/tests/test_integrator.py
@@ -63,7 +63,7 @@ class IntegratorTest(unittest.TestCase):
             0,
             R_SE * 2 * np.pi / T_E,
             0,
-            1 / M_E * M_E * R_SE * 2 * np.pi / T_E + 1 * R_L * 2 * np.pi / T_L,
+            R_SE * 2 * np.pi / T_E + 1 * R_L * 2 * np.pi / T_L,
         ])
 
         system = GravitationalSystem(r0=x0[:4],
@@ -78,11 +78,15 @@ class IntegratorTest(unittest.TestCase):
         ## plotting trajectories
         plt.figure()
         plt.plot(q[:, 0], q[:, 1], label="Earth")
+        plt.xlabel("x/AU")
+        plt.ylabel("y/AU")
         plt.legend()
         plt.show()
 
         plt.figure()
         plt.plot(q[:, 2] - q[:, 0], q[:, 3] - q[:, 1], label="Moon seen from Earth")
+        plt.xlabel("x/AU")
+        plt.ylabel("y/AU")
         plt.legend()
         plt.show()
 
@@ -102,8 +106,9 @@ class IntegratorTest(unittest.TestCase):
         self.assertTrue(np.greater(1e-10 + np.zeros(L.shape[0]), L - L[0]).all())
 
         ## checking error
-        dts = [dt, 2 * dt, 4 * dt]
-        errors = []
+        dts = [dt, 2 * dt, 4 * dt, 10 * dt, 50 * dt]
+        errors_E = []
+        errors_P = []
         ts = []
         for i in dts:
             system = GravitationalSystem(r0=x0[:4],
@@ -120,16 +125,35 @@ class IntegratorTest(unittest.TestCase):
             -G * M_S * M_E / np.linalg.norm(q_t[:,:2], axis=1) - G * M_S * M_L / np.linalg.norm(q_t[:,2:], axis=1) + \
             -G * M_E * M_L / np.linalg.norm(q_t[:,2:] - q_t[:,:2], axis=1)
 
-            errors.append((H - H[0]) / i**2)
+            P = p_t[:, :2] + p_t[:, 2:]
+            errors_E.append((H - H[0]) / i**2)
+            errors_P.append(np.linalg.norm(P, axis=1) - np.linalg.norm(P[0]))
             ts.append(t)
 
         plt.figure()
-        plt.plot(ts[0], errors[0], label="dt")
-        plt.plot(ts[1], errors[1], linestyle='--', label="2*dt")
-        plt.plot(ts[2], errors[2], linestyle=':', label="4*dt")
+        plt.plot(ts[0], errors_E[0], label="dt")
+        plt.plot(ts[1], errors_E[1], linestyle='--', label="2*dt")
+        plt.plot(ts[2], errors_E[2], linestyle=':', label="4*dt")
+        plt.plot(ts[3], errors_E[3], linestyle='-.', label="10*dt")
+        plt.plot(ts[4], errors_E[4], linestyle=':', label="50*dt")
         plt.xlabel("$t$")
         plt.ylabel("$\delta E(t)/(\Delta t)^2$")
         plt.legend()
+        plt.title("Energy error")
+        plt.tight_layout()
+        plt.show()
+
+        plt.figure()
+        plt.plot(ts[0], errors_P[0], label="dt")
+        plt.plot(ts[1], errors_P[1], linestyle='--', label="2*dt")
+        plt.plot(ts[2], errors_P[2], linestyle=':', label="4*dt")
+        plt.plot(ts[3], errors_P[3], linestyle='-.', label="10*dt")
+        plt.plot(ts[4], errors_P[4], linestyle=':', label="50*dt")
+        plt.xlabel("$t$")
+        plt.ylabel("$P(t)-P(0)$")
+        plt.legend()
+        plt.title("Linear momentum error")
+        plt.tight_layout()
         plt.show()
 
         ## checking time reversal: p -> -p
diff --git a/jobs/tests/test_jpl_data_query.py b/jobs/tests/test_jpl_data_query.py
index 5f884178f84e8577af62d260b53f1b646a074175..020ba0f60f0ce382edf9b7e1dac599e6497fffca 100644
--- a/jobs/tests/test_jpl_data_query.py
+++ b/jobs/tests/test_jpl_data_query.py
@@ -1,10 +1,5 @@
 """
 This unittest tests data query from JPL Horizons database 
-for the restricted three-body (sun-earth-moon) problem.
-The sun is fixed at the origin (center of mass). For 
-simplicity, the moon is assumed to be in the eliptic plane,
-so that vectors can be treated as two-dimensional.
-The trajectories are plotted in Heliographic Stonyhurst coordinates.
 """
 
 import unittest
@@ -21,22 +16,27 @@ time = parse_time('2024-01-01 00:00')
 class IntegratorTest(unittest.TestCase):
 
     def test_data_query(self):
-        moon = get_horizons_coord('301', {'start': time - 30 * u.day, 'stop': time + 30 * u.day, 'step': '1d'})
-
-        earth = get_horizons_coord('399', {'start': time - 30 * u.day, 'stop': time + 30 * u.day, 'step': '1d'})
-
-        sun = get_body_heliographic_stonyhurst('Sun', time)  # (lon, lat, radius) in (deg, deg, AU)
+        perihelion_14 = parse_time('2022-12-11 13:16')
+        psp = get_horizons_coord('Parker Solar Probe', {
+            'start': perihelion_14 - 50 * u.day,
+            'stop': perihelion_14 + 50 * u.day,
+            'step': '180m'
+        })
+        earth = get_body_heliographic_stonyhurst('Earth', perihelion_14)
 
         def coord_to_polar(coord):
             return coord.lon.to_value('rad'), coord.radius.to_value('AU')
 
         fig = plt.figure()
         ax = fig.add_subplot(projection='polar')
-        ax.plot(0, 0, 'o', label='Sun', color='orange', markersize=10)
-        ax.plot(*coord_to_polar(earth.transform_to(sun)), label='Earth', color='blue', linestyle='-.')
-
-        ax.plot(*coord_to_polar(moon.transform_to(sun)), label='Moon', color='purple', linestyle='dashed')
+        ax.plot(0, 0, 'o', label='Sun', color='orange')
+        ax.plot(*coord_to_polar(earth), 'o', label='Earth', color='blue')
+        ax.plot(*coord_to_polar(psp), label='PSP (as seen from Earth)', color='purple')
+        ax.plot(*coord_to_polar(psp.transform_to(earth)),
+                label='PSP (non-rotating frame)',
+                color='purple',
+                linestyle='dashed')
         ax.set_title('Stonyhurst heliographic coordinates')
-        ax.legend(loc='upper left', bbox_to_anchor=(1.0, 0.5))
+        ax.legend(loc='upper center')
 
         plt.show()
diff --git a/tasks/src/direct_simulation.py b/tasks/src/direct_simulation.py
index be29253e9e7f914d7bea61e9de0196549f6f7371..4d8b20a2bf09af7047614103c3be45ffbbb8f075 100644
--- a/tasks/src/direct_simulation.py
+++ b/tasks/src/direct_simulation.py
@@ -151,4 +151,4 @@ ax.set_xticks([])
 ax.set_yticks([])
 ax.set_zticks([])
 ax.legend()
-plt.show()
+plt.show()
\ No newline at end of file
diff --git a/tasks/src/stability_analysis.py b/tasks/src/stability_analysis.py
new file mode 100644
index 0000000000000000000000000000000000000000..b89082d1c6e4783ffddd226ee444dee7dffbaa7d
--- /dev/null
+++ b/tasks/src/stability_analysis.py
@@ -0,0 +1,204 @@
+"""
+This module performs stability analysis of the restricted 
+3-body system using realistic and synthetic initial data
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+
+from jobs.src.integrators import verlet
+from jobs.src.system import GravitationalSystem
+from tasks.src.utils import Scaler
+
+# system settings
+R_SE = 1  # distance between Earth and the Sun, 1 astronomical unit (AU)
+R_L = 2.57e-3 * R_SE  # distance between Earth and the Moon
+
+M_S = 1  # solar mass
+M_E = 3.00e-6 * M_S
+M_L = 3.69e-8 * M_S
+
+T_E = 1  # earth yr
+T_L = 27.3 / 365.3 * T_E
+
+G = 4 * np.pi**2  # AU^3 / m / yr^2
+
+# simulation settings
+T = 0.3  # yr
+n_order = 6
+dt = 4**(-n_order)  # yr
+
+x0 = np.array([
+    R_SE,
+    0,
+    R_SE + R_L,
+    0,
+    0,
+    R_SE * 2 * np.pi / T_E,
+    0,
+    R_SE * 2 * np.pi / T_E + 1 * R_L * 2 * np.pi / T_L,
+])
+SNR = [0.1, 0.05, 0.005]  # signal to noise ratio
+
+
+# force computation
+def force(q: np.ndarray) -> np.ndarray:
+    r1 = q[0:2]
+    r2 = q[2:4]
+
+    return np.array([
+        -G * M_E * (M_S * r1 / np.linalg.norm(r1)**3 + M_L * (r1 - r2) / np.linalg.norm(r1 - r2)**3),
+        -G * M_L * (M_S * r2 / np.linalg.norm(r2)**3 + M_E * (r2 - r1) / np.linalg.norm(r2 - r1)**3),
+    ]).flatten()
+
+
+system = GravitationalSystem(r0=x0[:4],
+                             v0=x0[4:],
+                             m=np.array([M_E, M_E, M_L, M_L]),
+                             t=np.linspace(0, T, int(T // dt)),
+                             force=force,
+                             solver=verlet)
+
+t, p_org, q_org = system.simulation()
+
+## Case 1: synthetic initial data
+# a) Gaussian noise added to position
+scaler_1_0 = Scaler(x0[:4], new_max=1.0, new_min=0.0)
+pos0_scaled = scaler_1_0.scale(x0[:4])
+
+for j in SNR:
+    q_pos_perturbed = []
+    for i in range(100):
+        pos0_perturbed = pos0_scaled + j * np.random.normal(size=4)
+        pos0_perturbed = scaler_1_0.inverse_scale(pos0_perturbed)
+        system = GravitationalSystem(r0=pos0_perturbed,
+                                     v0=x0[4:],
+                                     m=np.array([M_E, M_E, M_L, M_L]),
+                                     t=np.linspace(0, T, int(T // dt)),
+                                     force=force,
+                                     solver=verlet)
+        t, p, q = system.simulation()
+        q_pos_perturbed.append(q)
+
+    q_pos_perturbed = np.average(np.array(q_pos_perturbed), axis=0)
+
+    plt.figure()
+    plt.plot(q_org[:, 0], q_org[:, 1], label="Earth")
+    plt.plot(0, 0, 'o', label="Sun", color='pink', markersize=10)
+    plt.plot(q_pos_perturbed[:, 0], q_pos_perturbed[:, 1], label="av. perturbed Earth")
+    plt.xlabel("x/AU")
+    plt.ylabel("y/AU")
+    plt.grid()
+    plt.legend()
+    plt.title("Perturbed initial positions")
+    plt.tight_layout()
+
+    plt.savefig(f"pert_init_pos_earth_SNR_{j}.png")
+
+    plt.figure()
+    plt.plot(q_org[:, 2] - q_org[:, 0], q_org[:, 3] - q_org[:, 1], label="Moon")
+    plt.plot(0, 0, 'o', label="Sun", color='pink', markersize=10)
+    plt.plot(q_pos_perturbed[:, 2] - q_pos_perturbed[:, 0],
+             q_pos_perturbed[:, 3] - q_pos_perturbed[:, 1],
+             label="av. perturbed Moon")
+    plt.xlabel("x/AU")
+    plt.ylabel("y/AU")
+    plt.grid()
+    plt.legend()
+    plt.title("Perturbed initial positions")
+    plt.tight_layout()
+    plt.savefig(f"pert_init_pos_moon_SNR_{j}.png")
+
+# b) Gaussian noise added to velocity
+scaler_1_0 = Scaler(x0[4:], new_max=1.0, new_min=0.0)
+velo0_scaled = scaler_1_0.scale(x0[4:])
+
+for j in SNR:
+    q_velo_perturbed = []
+    for i in range(100):
+        velo0_perturbed = velo0_scaled + j * np.random.normal(size=4)
+        velo0_perturbed = scaler_1_0.inverse_scale(velo0_perturbed)
+        system = GravitationalSystem(r0=x0[:4],
+                                     v0=velo0_perturbed,
+                                     m=np.array([M_E, M_E, M_L, M_L]),
+                                     t=np.linspace(0, T, int(T // dt)),
+                                     force=force,
+                                     solver=verlet)
+        t, p, q = system.simulation()
+        q_velo_perturbed.append(q)
+
+    q_velo_perturbed = np.average(np.array(q_velo_perturbed), axis=0)
+
+    plt.figure()
+    plt.plot(q_org[:, 0], q_org[:, 1], label="Earth")
+    plt.plot(0, 0, 'o', label="Sun", color='pink', markersize=10)
+    plt.plot(q_velo_perturbed[:, 0], q_velo_perturbed[:, 1], label="av. perturbed Earth")
+    plt.xlabel("x/AU")
+    plt.ylabel("y/AU")
+    plt.grid()
+    plt.legend()
+    plt.title("Perturbed initial velocities")
+    plt.tight_layout()
+    plt.savefig(f"pert_init_velo_earth_SNR_{j}.png")
+
+    plt.figure()
+    plt.plot(q_org[:, 2] - q_org[:, 0], q_org[:, 3] - q_org[:, 1], label="Moon")
+    plt.plot(0, 0, 'o', label="Sun", color='pink', markersize=10)
+    plt.plot(q_velo_perturbed[:, 2] - q_velo_perturbed[:, 0],
+             q_velo_perturbed[:, 3] - q_velo_perturbed[:, 1],
+             label="av. perturbed Moon")
+    plt.xlabel("x/AU")
+    plt.ylabel("y/AU")
+    plt.grid()
+    plt.legend()
+    plt.title("Perturbed initial velocities")
+    plt.tight_layout()
+
+    plt.savefig(f"pert_init_velo_moon_SNR_{j}.png")
+
+# c) Gaussian noise added to position + velocity
+veloscaler_1_0 = Scaler(x0[4:], new_max=1.0, new_min=0.0)
+velo0_scaled = veloscaler_1_0.scale(x0[4:])
+posscaler_1_0 = Scaler(x0[:4], new_max=1.0, new_min=0.0)
+pos0_scaled = posscaler_1_0.scale(x0[:4])
+
+for j in SNR:
+    q_perturbed = []
+    for i in range(100):
+        pos0_perturbed = pos0_scaled + j * np.random.normal(size=4)
+        pos0_perturbed = posscaler_1_0.inverse_scale(pos0_perturbed)
+        velo0_perturbed = velo0_scaled + j * np.random.normal(size=4)
+        velo0_perturbed = veloscaler_1_0.inverse_scale(velo0_perturbed)
+        system = GravitationalSystem(r0=pos0_perturbed,
+                                     v0=velo0_perturbed,
+                                     m=np.array([M_E, M_E, M_L, M_L]),
+                                     t=np.linspace(0, T, int(T // dt)),
+                                     force=force,
+                                     solver=verlet)
+        t, p, q = system.simulation()
+        q_perturbed.append(q)
+
+    q_perturbed = np.average(np.array(q_perturbed), axis=0)
+
+    plt.figure()
+    plt.plot(q_org[:, 0], q_org[:, 1], label="Earth")
+    plt.plot(0, 0, 'o', label="Sun", color='pink', markersize=10)
+    plt.plot(q_perturbed[:, 0], q_perturbed[:, 1], label="av. perturbed Earth")
+    plt.xlabel("x/AU")
+    plt.ylabel("y/AU")
+    plt.grid()
+    plt.legend()
+    plt.title("Perturbed initial position + velocities")
+    plt.tight_layout()
+    plt.savefig(f"pert_init_pos_velo_earth_SNR_{j}.png")
+
+    plt.figure()
+    plt.plot(q_org[:, 2] - q_org[:, 0], q_org[:, 3] - q_org[:, 1], label="Moon")
+    plt.plot(0, 0, 'o', label="Sun", color='pink', markersize=10)
+    plt.plot(q_perturbed[:, 2] - q_perturbed[:, 0], q_perturbed[:, 3] - q_perturbed[:, 1], label="av. perturbed Moon")
+    plt.xlabel("x/AU")
+    plt.ylabel("y/AU")
+    plt.grid()
+    plt.legend()
+    plt.title("Perturbed initial position + velocities")
+    plt.tight_layout()
+    plt.savefig(f"pert_init_pos_velo_moon_SNR_{j}.png")
diff --git a/tasks/src/utils.py b/tasks/src/utils.py
index 6e3b6115a38dc0f95c73416c5b376ef68b897866..b87e134642c35ff8a42fd53743f2ba9b1d8ecd67 100644
--- a/tasks/src/utils.py
+++ b/tasks/src/utils.py
@@ -2,6 +2,22 @@ import numpy as np
 from astropy.coordinates import CartesianRepresentation, SphericalRepresentation
 
 
+class Scaler(object):
+    """
+    Scaler class for min-max scaling
+    """
+
+    def __init__(self, data, new_max, new_min):
+        self.own_min, self.own_max = np.min(data), np.max(data)
+        self.new_max, self.new_min = new_max, new_min
+
+    def scale(self, data: np.ndarray) -> np.ndarray:
+        return (data - self.own_min) / (self.own_max - self.own_min) * (self.new_max - self.new_min) + self.new_min
+
+    def inverse_scale(self, data: np.ndarray) -> np.ndarray:
+        return (data - self.new_min) * (self.own_max - self.own_min) / (self.new_max - self.new_min) + self.own_min
+
+
 def rotation_matrix(phi, theta) -> np.ndarray:
     """
     Rotate a vector represented in spherical coordinates.