|
| 1 | +import math |
| 2 | +import matplotlib.pyplot as plt |
| 3 | +from matplotlib.animation import FuncAnimation |
| 4 | + |
| 5 | +# Constants |
| 6 | +G = 6.674e-11 # Gravitational constant |
| 7 | +dt = 3600*24 # Time step (seconds) |
| 8 | + |
| 9 | +# Inputs |
| 10 | +mass_star = 1.989e30 # Mass of the star (e.g. the Sun) |
| 11 | +mass_planet = 5.972e24 # Mass of the planet (e.g. the Earth) |
| 12 | +num_steps = 2000 # Number of simulation steps |
| 13 | + |
| 14 | +# Initial Conditions |
| 15 | +planet_position = [1.5e11, 0] # Start at 1 AU (same as Earth) |
| 16 | +velocity = [3000, 25000] # Initial velocity (m/s) |
| 17 | +star_position = [0, 0] # Star is at the origin |
| 18 | + |
| 19 | +# Lists to store the path for plotting |
| 20 | +x_positions = [] |
| 21 | +y_positions = [] |
| 22 | +trail_positions = [] |
| 23 | +trail_lines = [] |
| 24 | + |
| 25 | +# Create the figure |
| 26 | +fig, ax = plt.subplots() |
| 27 | +planet_dot, = ax.plot([], [], 'bo') # Planet as blue dot |
| 28 | +star_dot, = ax.plot(0, 0, 'yo', label='Star') # Star as yellow dot |
| 29 | +orbit_path, = ax.plot([], [], 'b-', linewidth=0.5) # Planet's full path |
| 30 | +ax.set_xlim(-3e11, 3e11) |
| 31 | +ax.set_ylim(-3e11, 3e11) |
| 32 | +ax.set_aspect('equal') |
| 33 | +ax.grid(True) |
| 34 | +ax.set_title("Orbital Animation") |
| 35 | +ax.set_xlabel("X position (m)") |
| 36 | +ax.set_ylabel("Y position (m)") |
| 37 | +ax.legend() |
| 38 | + |
| 39 | +def update(frame): |
| 40 | + global planet_position, velocity, trail_lines |
| 41 | + |
| 42 | + # Compute distance components |
| 43 | + dx = star_position[0] - planet_position[0] |
| 44 | + dy = star_position[1] - planet_position[1] |
| 45 | + magnitude = math.sqrt(dx**2 + dy**2) |
| 46 | + |
| 47 | + # Compute gravitational force |
| 48 | + force = G * mass_star * mass_planet / magnitude**2 |
| 49 | + |
| 50 | + # Unit vector from planet to star |
| 51 | + unit_vector = [dx / magnitude, dy / magnitude] |
| 52 | + |
| 53 | + # Force vector |
| 54 | + force_vector = [force * unit_vector[0], force * unit_vector[1]] |
| 55 | + |
| 56 | + # Acceleration = F / m |
| 57 | + acceleration = [force_vector[0] / mass_planet, force_vector[1] / mass_planet] |
| 58 | + |
| 59 | + # Update velocity |
| 60 | + velocity[0] += acceleration[0] * dt |
| 61 | + velocity[1] += acceleration[1] * dt |
| 62 | + |
| 63 | + # Update position |
| 64 | + planet_position[0] += velocity[0] * dt |
| 65 | + planet_position[1] += velocity[1] * dt |
| 66 | + |
| 67 | + # Store position for full orbit |
| 68 | + x_positions.append(planet_position[0]) |
| 69 | + y_positions.append(planet_position[1]) |
| 70 | + |
| 71 | + # Store position for fading trail |
| 72 | + trail_positions.append(planet_position.copy()) |
| 73 | + max_trail_length = 100 |
| 74 | + if len(trail_positions) > max_trail_length: |
| 75 | + trail_positions.pop(0) |
| 76 | + |
| 77 | + # Remove previous trail lines |
| 78 | + for line in trail_lines: |
| 79 | + line.remove() |
| 80 | + trail_lines = [] |
| 81 | + |
| 82 | + # Draw fading trail |
| 83 | + for i in range(1, len(trail_positions)): |
| 84 | + alpha = i / len(trail_positions) |
| 85 | + line, = ax.plot( |
| 86 | + [trail_positions[i-1][0], trail_positions[i][0]], |
| 87 | + [trail_positions[i-1][1], trail_positions[i][1]], |
| 88 | + color='blue', |
| 89 | + alpha=alpha, |
| 90 | + linewidth=2 |
| 91 | + ) |
| 92 | + trail_lines.append(line) |
| 93 | + |
| 94 | + # Update planet and orbit |
| 95 | + planet_dot.set_data([planet_position[0]], [planet_position[1]]) |
| 96 | + orbit_path.set_data(x_positions, y_positions) |
| 97 | + |
| 98 | + return planet_dot, orbit_path, *trail_lines |
| 99 | + |
| 100 | +# Animate! |
| 101 | +ani = FuncAnimation(fig, update, frames=num_steps, interval=10, blit=False) |
| 102 | +plt.show() |
0 commit comments