Funny plots

.gitlab-ci.yml

@@ -17,18 +17,6 @@ test-jobs:
     - export PYTHONPATH="$PYTHONPATH:$CI_PROJECT_DIR/jobs/src"
     - poetry run pytest jobs/tests/
 
-test-tasks:
-  stage: test-tasks
-  script:
-    - pip install poetry
-    - poetry --version
-    - cd build/
-    - poetry install --no-root
-    - source $(poetry env info --path)/bin/activate
-    - cd ..
-    - export PYTHONPATH="$PYTHONPATH:$CI_PROJECT_DIR/tasks/src"
-    - poetry run pytest tasks/tests/
-  
 pages:
   stage: build-page
   script:

build/pyproject.toml

@@ -6,7 +6,7 @@ authors = []
 readme = "README.md"
 
 [tool.poetry.dependencies]
-python = "3.11.8"
+python = "~3.11"
 numpy = "1.24.4"
 matplotlib = "3.8.0"
 astropy = "6.0.0"

jobs/src/bht_algorithm_3D.py

@@ -7,7 +7,7 @@ import numpy as np
 
 class MainApp:
 
-    def __init__(self, particles):
+    def __init__(self):
         """Initialize the MainApp with a root node."""
         self.rootNode = TreeNode(x=0, y=0, z=0, width=200, height=200, depth=200)
 

pyproject.toml

@@ -6,7 +6,7 @@ authors = []
 readme = "README.md"
 
 [tool.poetry.dependencies]
-python = "3.11.8"
+python = "~3.11"
 toml = "^0.10.2"
 
 [tool.poetry.group.dev.dependencies]

tasks/src/direct_simulation.py

@@ -10,6 +10,7 @@ import numpy as np
 from sunpy.coordinates import (get_body_heliographic_stonyhurst, get_horizons_coord)
 from sunpy.time import parse_time
 
+from jobs.src.bht_algorithm_3D import MainApp, Particle
 from jobs.src.integrators import verlet
 from jobs.src.system import GravitationalSystem
 from tasks.src.utils import (coord_to_polar, xyz_to_polar)
@@ -77,15 +78,51 @@ system = GravitationalSystem(r0=x0[:6],
 
 t, p, q = system.simulation()
 
+
+# run bht simulation
+# force computation
+def force(q: np.ndarray) -> np.ndarray:
+    r1 = q[0:3]  # Earth coordinates
+    r2 = q[3:6]  # Moon coordinates
+
+    sun = Particle(0, 0, 0, M_S)
+    earth = Particle(r1[0], r1[1], r1[2], M_E)
+    moon = Particle(r2[0], r2[1], r2[2], M_L)
+
+    particles = [sun, earth, moon]
+
+    barnes_hut = MainApp()  # Initialize Barnes-Hut algorithm instance
+    barnes_hut.BuildTree(particles)  # Build the Barnes-Hut tree with particles
+    barnes_hut.rootNode.ComputeMassDistribution()  #Compute the center of mass of the tree nodes
+
+    f_earth = barnes_hut.rootNode.CalculateForceFromTree(earth)
+    f_moon = barnes_hut.rootNode.CalculateForceFromTree(moon)
+
+    return np.concatenate((f_earth, f_moon), axis=0)
+
+
+system = GravitationalSystem(r0=x0[:6],
+                             v0=x0[6:],
+                             m=np.array([M_E, M_E, M_E, M_L, M_L, M_L]),
+                             t=np.linspace(0, T, int(T // dt)),
+                             force=force,
+                             solver=verlet)
+
+t, p_bht, q_bht = system.simulation()
+
 # plot trajectory in polar coordinates
 q_E_polar = np.array([xyz_to_polar(pos[0], pos[1], pos[2]) for pos in q])
 q_M_polar = np.array([xyz_to_polar(pos[3], pos[4], pos[5]) for pos in q])
+q_E_bht_polar = np.array([xyz_to_polar(pos[0], pos[1], pos[2]) for pos in q_bht])
+q_M_bht_polar = np.array([xyz_to_polar(pos[3], pos[4], pos[5]) for pos in q_bht])
 
 fig = plt.figure()
 ax = fig.add_subplot(projection='polar')
-ax.plot(0, 0, 'o', label='Sun', color='orange', markersize=10)
+ax.plot(0, 0, 'o', label='Sun', color='pink', markersize=10)
 ax.plot(q_E_polar[:, 0], q_E_polar[:, 1], label='Earth', color='blue', linestyle='-.')
 ax.plot(q_M_polar[:, 0], q_M_polar[:, 1], label='Moon', color='purple', linestyle='dashed')
+ax.plot(q_E_bht_polar[:, 0], q_E_bht_polar[:, 1], label='Earth(BHT)', linestyle='-.')
+ax.plot(q_M_bht_polar[:, 0], q_M_bht_polar[:, 1], label='Moon(BHT)', linestyle='dashed')
 ax.plot(*coord_to_polar(earth.transform_to(sun)), label='Earth(JPL)', color='green', linestyle='-.')
 ax.plot(*coord_to_polar(moon.transform_to(sun)), label='Moon(JPL)', color='red', linestyle='dashed')
 ax.set_title('Stonyhurst heliographic coordinates')
@@ -100,9 +137,11 @@ earth_JPL = np.array([(i.x.value, i.y.value, i.z.value) for i in earth_xyz])
 
 fig = plt.figure()
 ax = fig.add_subplot(projection='3d')
-ax.plot(0, 0, 0, 'o', label='Sun', color='orange', markersize=10)
+ax.plot(0, 0, 0, 'o', label='Sun', color='pink', markersize=10)
 ax.plot(q[:, 0], q[:, 1], q[:, 2], label="Earth")
 ax.plot(q[:, 3], q[:, 4], q[:, 5], label="Moon")
+ax.plot(q_bht[:, 0], q_bht[:, 1], q_bht[:, 2], label="Earth(BHT)")
+ax.plot(q_bht[:, 3], q_bht[:, 4], q_bht[:, 5], label="Moon(BHT)")
 ax.plot(earth_JPL[:, 0], earth_JPL[:, 1], earth_JPL[:, 2], '--', label='Earth(JPL)')
 ax.plot(moon_JPL[:, 0], moon_JPL[:, 1], moon_JPL[:, 2], '-.', label='Moon(JPL)')
 ax.set_xlabel("X")