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()