Skip to content

Commit 4590147

Browse files
New project!
1 parent 3a9a237 commit 4590147

1 file changed

Lines changed: 112 additions & 0 deletions

File tree

  • multi-body-orbital-mechanics
Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
import matplotlib.pyplot as plt
2+
import matplotlib.animation as animation
3+
import math
4+
5+
# --- Physics constants ---
6+
G = 6.67430e-11 # gravitational constant
7+
dt = 60 * 60 * 6 # time step (6 hours)
8+
frame_count = 0 # to animate glow pulsing
9+
10+
# --- Planet data ---
11+
planets = [
12+
{
13+
"pos": [0.39 * 1.496e11, 0], # Mercury
14+
"vel": [0, 47360],
15+
"mass": 3.285e23,
16+
"x_path": [],
17+
"y_path": [],
18+
"color": "gray"
19+
},
20+
{
21+
"pos": [0.72 * 1.496e11, 0], # Venus
22+
"vel": [0, 35020],
23+
"mass": 4.867e24,
24+
"x_path": [],
25+
"y_path": [],
26+
"color": "orange"
27+
},
28+
{
29+
"pos": [1.496e11, 0], # Earth
30+
"vel": [0, 29780],
31+
"mass": 5.972e24,
32+
"x_path": [],
33+
"y_path": [],
34+
"color": "blue"
35+
},
36+
{
37+
"pos": [1.52 * 1.496e11, 0], # Mars
38+
"vel": [0, 24077],
39+
"mass": 6.417e23,
40+
"x_path": [],
41+
"y_path": [],
42+
"color": "red"
43+
}
44+
]
45+
46+
# --- Setup figure ---
47+
fig, ax = plt.subplots()
48+
ax.set_facecolor("black")
49+
ax.set_aspect('equal', adjustable='box')
50+
51+
def update_positions():
52+
for i, p1 in enumerate(planets):
53+
fx = fy = 0
54+
for j, p2 in enumerate(planets):
55+
if i != j:
56+
dx = p2["pos"][0] - p1["pos"][0]
57+
dy = p2["pos"][1] - p1["pos"][1]
58+
dist = math.sqrt(dx**2 + dy**2)
59+
F = G * p1["mass"] * p2["mass"] / dist**2
60+
fx += F * dx / dist
61+
fy += F * dy / dist
62+
63+
# Add Sun gravity
64+
dx = 0 - p1["pos"][0]
65+
dy = 0 - p1["pos"][1]
66+
dist = math.sqrt(dx**2 + dy**2)
67+
F = G * p1["mass"] * 1.989e30 / dist**2
68+
fx += F * dx / dist
69+
fy += F * dy / dist
70+
71+
# Update velocity & position
72+
p1["vel"][0] += fx / p1["mass"] * dt
73+
p1["vel"][1] += fy / p1["mass"] * dt
74+
p1["pos"][0] += p1["vel"][0] * dt
75+
p1["pos"][1] += p1["vel"][1] * dt
76+
77+
# Save path
78+
p1["x_path"].append(p1["pos"][0])
79+
p1["y_path"].append(p1["pos"][1])
80+
81+
def animate(frame):
82+
global frame_count
83+
frame_count += 1
84+
ax.clear()
85+
ax.set_facecolor("black")
86+
ax.set_aspect('equal', adjustable='box')
87+
ax.set_xlim(-2.5e11, 2.5e11)
88+
ax.set_ylim(-2.5e11, 2.5e11)
89+
ax.axis('off')
90+
91+
update_positions()
92+
93+
# --- Sun breathing glow ---
94+
glow_alpha = 0.3 + 0.25 * (math.sin(frame_count * 0.05) + 1) / 2
95+
ax.scatter(0, 0, color='yellow', s=200, zorder=3)
96+
ax.scatter(0, 0, color='yellow', s=900, alpha=glow_alpha, zorder=1)
97+
98+
# --- Planets ---
99+
for p in planets:
100+
# Orbit path
101+
ax.plot(p["x_path"], p["y_path"], color=p["color"], lw=1)
102+
103+
# Planet
104+
ax.scatter(p["pos"][0], p["pos"][1], color=p["color"], s=30, zorder=3)
105+
106+
# Glow halo
107+
planet_alpha = 0.3 + 0.25 * (math.sin(frame_count * 0.1) + 1) / 2
108+
ax.scatter(p["pos"][0], p["pos"][1],
109+
color=p["color"], s=200, alpha=planet_alpha, zorder=1)
110+
111+
ani = animation.FuncAnimation(fig, animate, frames=360, interval=20)
112+
plt.show()

0 commit comments

Comments
 (0)