Skip to content
Snippets Groups Projects
Commit 5ee7d1a3 authored by nguyed99's avatar nguyed99
Browse files

Funny plots

parent 88caa309
Branches
No related tags found
No related merge requests found
Pipeline #59628 passed
......@@ -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:
......
......@@ -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"
......
......@@ -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)
......
......@@ -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]
......
......@@ -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")
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment