51 lines
1.1 KiB
Python
51 lines
1.1 KiB
Python
import numpy as np
|
|
|
|
import matplotlib.pyplot as plt
|
|
|
|
|
|
|
|
# Constants
|
|
G = 6.67430e-11 # Gravitational constant (m^3 kg^-1 s^-2)
|
|
M = 5.972e24 # Mass of Earth (kg)
|
|
R = 6.371e6 # Radius of Earth (m)
|
|
|
|
# Initial conditions
|
|
r0 = R + 400e3 # Initial distance from Earth's center (400 km altitude)
|
|
v0 = np.sqrt(G * M / r0) # Circular orbit velocity
|
|
|
|
# Simulation parameters
|
|
dt = 10 # Time step (s)
|
|
num_steps = 10000
|
|
|
|
# Arrays to store position
|
|
x = np.zeros(num_steps)
|
|
y = np.zeros(num_steps)
|
|
vx = np.zeros(num_steps)
|
|
vy = np.zeros(num_steps)
|
|
|
|
# Initial position and velocity
|
|
x[0] = r0
|
|
y[0] = 0
|
|
vx[0] = 0
|
|
vy[0] = v0
|
|
|
|
for i in range(1, num_steps):
|
|
r = np.sqrt(x[i-1]**2 + y[i-1]**2)
|
|
ax = -G * M * x[i-1] / r**3
|
|
ay = -G * M * y[i-1] / r**3
|
|
vx[i] = vx[i-1] + ax * dt
|
|
vy[i] = vy[i-1] + ay * dt
|
|
x[i] = x[i-1] + vx[i] * dt
|
|
y[i] = y[i-1] + vy[i] * dt
|
|
|
|
# Plotting
|
|
plt.figure(figsize=(6,6))
|
|
earth = plt.Circle((0, 0), R, color='b', alpha=0.3)
|
|
plt.gca().add_patch(earth)
|
|
plt.plot(x, y, label='Orbit')
|
|
plt.xlabel('x (m)')
|
|
plt.ylabel('y (m)')
|
|
plt.title('Orbital Mechanics Demonstration')
|
|
plt.axis('equal')
|
|
plt.legend()
|
|
plt.show() |