The movements of planets in our solar system have fascinated humanity for millennia. Today, we can use Python to calculate and visualise the positions of planets with impressive accuracy. In this article, we'll build a program that calculates the positions of planets for any given date and renders a beautiful vector graphic representation of the solar system.
Our approach will combine astronomical calculations with the powerful visualisation capabilities of Python libraries. We'll learn about astronomical coordinate systems, orbital mechanics basics, and how to transform complex calculations into an elegant visual model.
Understanding Planetary Motion
Planets move around the Sun in elliptical orbits. The exact position of a planet at any given time can be determined using Kepler's laws of planetary motion and a set of orbital elements. While the complete calculations involve complex mathematics, we can use simplified models that provide good approximations for visualisation purposes.
For each planet, we'll need the following orbital elements:
- Semi-major axis (a): The average distance from the planet to the Sun
- Eccentricity (e): How elliptical (rather than circular) the orbit is
- Inclination (i): The tilt of the orbital plane relative to the Earth's orbital plane
- Longitude of the ascending node (Ω): Where the orbit crosses the reference plane
- Argument of perihelion (ω): The orientation of the ellipse in the orbital plane
- Mean anomaly (M₀): The position of the planet at a reference time
Building Our Solar System Model
Let's start by setting up our project and importing the necessary libraries. We'll use matplotlib for visualisation, numpy for numerical calculations, and datetime for handling dates.
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import numpy as np
from datetime import datetime, timedelta
import math
Next, we'll define the orbital elements for each planet. These values are approximations that change slightly over time, but they're accurate enough for our visualisation.
# Orbital elements of planets (Mercury to Neptune)
# Format: [semi-major axis (AU), eccentricity, inclination (deg),
# longitude of ascending node (deg), argument of perihelion (deg),
# mean anomaly at J2000 (deg), orbital period (years)]
PLANETS = {
"Mercury": [0.387098, 0.205630, 7.005, 48.331, 29.124, 174.795, 0.240846],
"Venus": [0.723332, 0.006772, 3.395, 76.680, 54.884, 50.416, 0.615197],
"Earth": [1.000000, 0.016709, 0.000, -11.260, 114.208, 358.617, 1.000174],
"Mars": [1.523679, 0.093396, 1.851, 49.558, 286.502, 19.373, 1.880765],
"Jupiter": [5.204267, 0.048775, 1.305, 100.492, 273.867, 20.020, 11.862615],
"Saturn": [9.582069, 0.055723, 2.484, 113.642, 339.392, 317.020, 29.655475],
"Uranus": [19.229412, 0.046321, 0.770, 74.000, 96.998, 142.238, 84.039492],
"Neptune": [30.103658, 0.009456, 1.769, 131.784, 276.336, 256.228, 164.793623]
}
# Planet colours for visualisation
PLANET_COLORS = {
"Mercury": "#8C8680",
"Venus": "#E8BB67",
"Earth": "#4F6bed",
"Mars": "#CD5C5C",
"Jupiter": "#E0CFAF",
"Saturn": "#F0E6AA",
"Uranus": "#9AECF9",
"Neptune": "#4169E1"
}
# Relative size of planets (for visualisation, not to scale)
PLANET_SIZES = {
"Mercury": 0.38,
"Venus": 0.95,
"Earth": 1.0,
"Mars": 0.53,
"Jupiter": 11.2,
"Saturn": 9.45,
"Uranus": 4.0,
"Neptune": 3.88
}
Calculating Planetary Positions
Now, let's write the functions to calculate the position of a planet at a given date. We'll convert orbital elements to heliocentric (Sun-centered) Cartesian coordinates.
def calculate_planet_position(planet_name, date):
"""Calculate the position of a planet on a given date."""
# Get orbital elements
a, e, i, node, peri, M0, period = PLANETS[planet_name]
# Convert angles from degrees to radians
i = math.radians(i)
node = math.radians(node)
peri = math.radians(peri)
# Calculate days since J2000 (January 1, 2000 12:00 UTC)
j2000 = datetime(2000, 1, 1, 12, 0, 0)
days_since_j2000 = (date - j2000).total_seconds() / (24 * 3600)
# Calculate the mean anomaly at the given date
# M0 is in degrees, need to convert
M0_rad = math.radians(M0)
n = 2 * math.pi / (period * 365.25) # Mean motion (radians per day)
M = M0_rad + n * days_since_j2000
# Solve Kepler's equation using Newton-Raphson method
E = M # Initial guess for eccentric anomaly
for _ in range(5): # Usually converges quickly
E = E - (E - e * math.sin(E) - M) / (1 - e * math.cos(E))
# Calculate true anomaly
v = 2 * math.atan2(math.sqrt(1 + e) * math.sin(E/2), math.sqrt(1 - e) * math.cos(E/2))
# Calculate heliocentric distance
r = a * (1 - e * math.cos(E))
# Calculate heliocentric Cartesian coordinates
# First, in the orbital plane
x_orb = r * math.cos(v)
y_orb = r * math.sin(v)
# Then, apply rotations to get coordinates in the reference plane
x = (x_orb * (math.cos(peri) * math.cos(node) - math.sin(peri) * math.sin(node) * math.cos(i)) -
y_orb * (math.sin(peri) * math.cos(node) + math.cos(peri) * math.sin(node) * math.cos(i)))
y = (x_orb * (math.cos(peri) * math.sin(node) + math.sin(peri) * math.cos(node) * math.cos(i)) +
y_orb * (math.cos(peri) * math.cos(node) * math.cos(i) - math.sin(peri) * math.sin(node)))
z = (x_orb * math.sin(peri) * math.sin(i) + y_orb * math.cos(peri) * math.sin(i))
return x, y, z
This function implements Kepler's equations to calculate planetary positions. Here's how it works:
- First, we convert all angles from degrees to radians for mathematical calculations.
- We calculate the number of days that have passed since J2000 (a standard astronomical epoch).
- The mean anomaly (M) is calculated by adding the planet's motion since J2000 to its initial position.
- We solve Kepler's equation iteratively to find the eccentric anomaly (E).
- From E, we calculate the true anomaly (v) and the distance from the Sun (r).
- Finally, we apply rotation matrices to transform orbital plane coordinates to the reference plane coordinates.
Visualising the Solar System
Now that we can calculate planetary positions, let's create functions to visualise the solar system using matplotlib:
def plot_solar_system(date, view="top", max_au=32, show_names=True):
"""
Plot the solar system for a given date.
Parameters:
- date: datetime object
- view: 'top' for top-down view, 'side' for side view
- max_au: maximum distance from Sun to show (in AU)
- show_names: whether to show planet names
"""
fig, ax = plt.subplots(figsize=(12, 12))
# Set plot limits
ax.set_xlim(-max_au, max_au)
ax.set_ylim(-max_au, max_au)
# Make plot square
ax.set_aspect('equal')
# Draw the Sun
sun = Circle((0, 0), 0.25, color='yellow', zorder=10)
ax.add_artist(sun)
# For each planet
for planet in PLANETS:
# Get position
x, y, z = calculate_planet_position(planet, date)
# Choose coordinates based on view
if view == "top":
plot_x, plot_y = x, y
elif view == "side":
plot_x, plot_y = x, z
# Draw orbit - use actual semi-major axis for correct scaling
a = PLANETS[planet][0] # Semi-major axis in AU
theta = np.linspace(0, 2 * np.pi, 100)
# Create proper elliptical orbit
e = PLANETS[planet][1] # Eccentricity
if view == "top":
# For top view, we can show proper elliptical orbits
orbit_x = a * np.cos(theta)
orbit_y = a * np.sin(theta)
# Apply eccentricity (simplified)
center_x = -a * e
orbit_x = orbit_x + center_x
else:
# For side view, we'll use circular approximation since
# we're not accounting for all 3D rotations
orbit_x = a * np.cos(theta)
orbit_y = a * np.sin(theta) * math.sin(math.radians(PLANETS[planet][2]))
# Draw orbit path (with transparency)
ax.plot(orbit_x, orbit_y, color='gray', alpha=0.3)
# Draw actual planet
size = 0.2 * max(0.3, min(1.0, np.log10(PLANET_SIZES[planet]) + 0.5))
planet_circle = Circle((plot_x, plot_y), size,
color=PLANET_COLORS[planet], zorder=5)
ax.add_artist(planet_circle)
# Label planets if requested
if show_names:
ax.annotate(planet, (plot_x, plot_y), xytext=(5, 5),
textcoords='offset points', fontsize=10, color='white')
# Add distance indicators (circular grid lines)
for distance in [5, 10, 15, 20, 25, 30]:
if distance <= max_au:
circle = plt.Circle((0, 0), distance, color="gray", fill=False, alpha=0.2, linestyle="--")
ax.add_artist(circle)
# Label the distance grid
ax.annotate(f"{distance} AU", (distance, 0), color="gray", alpha=0.7,
fontsize=8, ha='left', va='bottom')
# Set title and labels
date_str = date.strftime("%Y-%m-%d")
ax.set_title(f"Solar System on {date_str} ({view}-view)", fontsize=16, color='white')
ax.set_xlabel("x (AU)" if view == "top" else "x (AU)", color='white')
ax.set_ylabel("y (AU)" if view == "top" else "z (AU)", color='white')
# Grid
ax.grid(True, alpha=0.3)
# Remove ticks for cleaner look
ax.tick_params(axis='both', which='both', length=0, colors='white')
# Background
ax.set_facecolor('black')
return fig
This visualisation function handles a few important tasks:
- Sets up a square plot with the Sun at the center
- Calculates and plots the position of each planet
- Draws approximate orbital paths
- Provides both top-down and side views of the solar system
- Labels planets and adds a title with the date
While this is a simplified visualisation (real orbits are elliptical and tilted), it still provides an accurate representation of the relative positions of planets at any given date.
Interactive Exploration
Let's add an interactive function that allows us to visualise planetary positions on any date:
def explore_solar_system(date_str=None, days_from_now=0):
"""
Explore the solar system on a specific date or days from today.
Parameters:
- date_str: Date string in 'YYYY-MM-DD' format
- days_from_now: Number of days from today (used if date_str is None)
"""
if date_str:
try:
date = datetime.strptime(date_str, "%Y-%m-%d")
except ValueError:
print("Invalid date format. Please use YYYY-MM-DD.")
return
else:
date = datetime.now() + timedelta(days=days_from_now)
# Create top and side views
fig_top = plot_solar_system(date, view="top")
fig_side = plot_solar_system(date, view="side")
plt.tight_layout()
plt.show()
# Print planetary positions
print(f"Planetary positions on {date.strftime('%Y-%m-%d')}:")
print("-" * 50)
print(f"{'Planet':<10} {'X (AU)':<10} {'Y (AU)':<10} {'Z (AU)':<10}")
print("-" * 50)
for planet in PLANETS:
x, y, z = calculate_planet_position(planet, date)
print(f"{planet:<10} {x:<10.3f} {y:<10.3f} {z:<10.3f}")
# Example usage:
# explore_solar_system(date_str="2023-12-25") # Christmas 2023
# explore_solar_system(days_from_now=365) # One year from now
With this interactive function, we can now explore the solar system on any date we choose, seeing both top-down and side views, plus the exact coordinates of each planet.
Adding Animation
To make our visualisation even more engaging, let's create an animation that shows the planets orbiting the Sun over time:
import matplotlib.animation as animation
def animate_solar_system(start_date, days=365, interval=20):
"""
Create an animation of the solar system over a period of time.
Parameters:
- start_date: datetime object for the starting date
- days: number of days to animate
- interval: time between frames in milliseconds
"""
# Set up the figure
fig, ax = plt.subplots(figsize=(10, 10))
max_au = 32 # Maximum distance to show
# Function to initialize the plot
def init():
ax.clear()
ax.set_xlim(-max_au, max_au)
ax.set_ylim(-max_au, max_au)
ax.set_aspect('equal')
ax.set_facecolor('black')
return []
# Function to update the plot for each frame
def update(frame):
ax.clear()
ax.set_xlim(-max_au, max_au)
ax.set_ylim(-max_au, max_au)
ax.set_aspect('equal')
ax.set_facecolor('black')
# Current date
current_date = start_date + timedelta(days=frame)
date_str = current_date.strftime("%Y-%m-%d")
ax.set_title(f"Solar System on {date_str}", fontsize=16, color='white')
# Draw the Sun
sun = Circle((0, 0), 0.25, color='yellow', zorder=10)
ax.add_artist(sun)
# Draw planets and their orbits
objects = []
for planet in PLANETS:
# Get position
x, y, z = calculate_planet_position(planet, current_date)
# Draw orbit (simplified circular orbit)
a = PLANETS[planet][0]
theta = np.linspace(0, 2 * np.pi, 100)
orbit_x = a * np.cos(theta)
orbit_y = a * np.sin(theta)
ax.plot(orbit_x, orbit_y, color='gray', alpha=0.3)
# Draw planet
size = 0.2 * max(0.3, min(1.0, np.log10(PLANET_SIZES[planet]) + 0.5))
planet_circle = Circle((x, y), size, color=PLANET_COLORS[planet], zorder=5)
ax.add_artist(planet_circle)
objects.append(planet_circle)
# Label planet
ax.annotate(planet, (x, y), xytext=(5, 5),
textcoords='offset points', fontsize=10, color='white')
ax.grid(True, alpha=0.3)
ax.set_xlabel("x (AU)", color='white')
ax.set_ylabel("y (AU)", color='white')
ax.tick_params(axis='both', colors='white')
return objects
# Create the animation
ani = animation.FuncAnimation(fig, update, frames=days,
init_func=init, blit=False, interval=interval)
plt.tight_layout()
return ani
# Example usage:
# start_date = datetime(2023, 1, 1)
# ani = animate_solar_system(start_date, days=365*2) # Animate 2 years
# ani.save('solar_system.mp4', writer='ffmpeg', fps=30, dpi=100)
# plt.show()
This animation function creates a video showing the movement of planets around the Sun over a specified time period. You can save the animation to a file or display it interactively.
Putting It All Together
Let's create a main function that brings everything together and provides a complete interface to our solar system visualisation:
def main():
"""Main function to run the solar system visualisation program."""
print("Solar System visualisation")
print("=" * 30)
while True:
print("\nWhat would you like to do?")
print("1. View solar system on a specific date")
print("2. View solar system today")
print("3. Create animation")
print("4. Exit")
choice = input("\nEnter your choice (1-4): ")
if choice == '1':
date_str = input("Enter date (YYYY-MM-DD): ")
explore_solar_system(date_str=date_str)
elif choice == '2':
explore_solar_system()
elif choice == '3':
start_date_str = input("Enter start date (YYYY-MM-DD, default=today): ")
if start_date_str:
try:
start_date = datetime.strptime(start_date_str, "%Y-%m-%d")
except ValueError:
print("Invalid date format. Using today's date.")
start_date = datetime.now()
else:
start_date = datetime.now()
days_str = input("Enter number of days to animate (default=365): ")
days = int(days_str) if days_str else 365
print("Creating animation (this may take a while)...")
ani = animate_solar_system(start_date, days=days)
save = input("Save animation to file? (y/n, default=n): ")
if save.lower() == 'y':
filename = input("Enter filename (default=solar_system.mp4): ")
if not filename:
filename = "solar_system.mp4"
ani.save(filename, writer='ffmpeg', fps=30, dpi=100)
print(f"Animation saved to {filename}")
plt.show()
elif choice == '4':
break
else:
print("Invalid choice. Please try again.")
if __name__ == "__main__":
main()
The Complete Code
Here's the complete program that brings together all the components we've discussed:
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import numpy as np
from datetime import datetime, timedelta
import math
import matplotlib.animation as animation
# Orbital elements of planets (Mercury to Neptune)
# Format: [semi-major axis (AU), eccentricity, inclination (deg),
# longitude of ascending node (deg), argument of perihelion (deg),
# mean anomaly at J2000 (deg), orbital period (years)]
PLANETS = {
"Mercury": [0.387098, 0.205630, 7.005, 48.331, 29.124, 174.795, 0.240846],
"Venus": [0.723332, 0.006772, 3.395, 76.680, 54.884, 50.416, 0.615197],
"Earth": [1.000000, 0.016709, 0.000, -11.260, 114.208, 358.617, 1.000174],
"Mars": [1.523679, 0.093396, 1.851, 49.558, 286.502, 19.373, 1.880765],
"Jupiter": [5.204267, 0.048775, 1.305, 100.492, 273.867, 20.020, 11.862615],
"Saturn": [9.582069, 0.055723, 2.484, 113.642, 339.392, 317.020, 29.655475],
"Uranus": [19.229412, 0.046321, 0.770, 74.000, 96.998, 142.238, 84.039492],
"Neptune": [30.103658, 0.009456, 1.769, 131.784, 276.336, 256.228, 164.793623]
}
# Planet colours for visualisation
PLANET_COLORS = {
"Mercury": "#8C8680",
"Venus": "#E8BB67",
"Earth": "#4F6bed",
"Mars": "#CD5C5C",
"Jupiter": "#E0CFAF",
"Saturn": "#F0E6AA",
"Uranus": "#9AECF9",
"Neptune": "#4169E1"
}
# Relative size of planets (for visualisation, not to scale)
PLANET_SIZES = {
"Mercury": 0.38,
"Venus": 0.95,
"Earth": 1.0,
"Mars": 0.53,
"Jupiter": 11.2,
"Saturn": 9.45,
"Uranus": 4.0,
"Neptune": 3.88
}
def calculate_planet_position(planet_name, date):
"""Calculate the position of a planet on a given date."""
# Get orbital elements
a, e, i, node, peri, M0, period = PLANETS[planet_name]
# Convert angles from degrees to radians
i = math.radians(i)
node = math.radians(node)
peri = math.radians(peri)
# Calculate days since J2000 (January 1, 2000 12:00 UTC)
j2000 = datetime(2000, 1, 1, 12, 0, 0)
days_since_j2000 = (date - j2000).total_seconds() / (24 * 3600)
# Calculate the mean anomaly at the given date
# M0 is in degrees, need to convert
M0_rad = math.radians(M0)
n = 2 * math.pi / (period * 365.25) # Mean motion (radians per day)
M = M0_rad + n * days_since_j2000
# Solve Kepler's equation using Newton-Raphson method
E = M # Initial guess for eccentric anomaly
for _ in range(5): # Usually converges quickly
E = E - (E - e * math.sin(E) - M) / (1 - e * math.cos(E))
# Calculate true anomaly
v = 2 * math.atan2(math.sqrt(1 + e) * math.sin(E/2), math.sqrt(1 - e) * math.cos(E/2))
# Calculate heliocentric distance
r = a * (1 - e * math.cos(E))
# Calculate heliocentric Cartesian coordinates
# First, in the orbital plane
x_orb = r * math.cos(v)
y_orb = r * math.sin(v)
# Then, apply rotations to get coordinates in the reference plane
x = (x_orb * (math.cos(peri) * math.cos(node) - math.sin(peri) * math.sin(node) * math.cos(i)) -
y_orb * (math.sin(peri) * math.cos(node) + math.cos(peri) * math.sin(node) * math.cos(i)))
y = (x_orb * (math.cos(peri) * math.sin(node) + math.sin(peri) * math.cos(node) * math.cos(i)) +
y_orb * (math.cos(peri) * math.cos(node) * math.cos(i) - math.sin(peri) * math.sin(node)))
z = (x_orb * math.sin(peri) * math.sin(i) + y_orb * math.cos(peri) * math.sin(i))
return x, y, z
def plot_solar_system(date, view="top", max_au=32, show_names=True):
"""
Plot the solar system for a given date.
Parameters:
- date: datetime object
- view: 'top' for top-down view, 'side' for side view
- max_au: maximum distance from Sun to show (in AU)
- show_names: whether to show planet names
"""
fig, ax = plt.subplots(figsize=(12, 12))
# Set plot limits
ax.set_xlim(-max_au, max_au)
ax.set_ylim(-max_au, max_au)
# Make plot square
ax.set_aspect('equal')
# Draw the Sun
sun = Circle((0, 0), 0.25, color='yellow', zorder=10)
ax.add_artist(sun)
# For each planet
for planet in PLANETS:
# Get position
x, y, z = calculate_planet_position(planet, date)
# Choose coordinates based on view
if view == "top":
plot_x, plot_y = x, y
elif view == "side":
plot_x, plot_y = x, z
# Draw orbit - use actual semi-major axis for correct scaling
a = PLANETS[planet][0] # Semi-major axis in AU
theta = np.linspace(0, 2 * np.pi, 100)
# Create proper elliptical orbit
e = PLANETS[planet][1] # Eccentricity
if view == "top":
# For top view, we can show proper elliptical orbits
orbit_x = a * np.cos(theta)
orbit_y = a * np.sin(theta)
# Apply eccentricity (simplified)
center_x = -a * e
orbit_x = orbit_x + center_x
else:
# For side view, we'll use circular approximation since
# we're not accounting for all 3D rotations
orbit_x = a * np.cos(theta)
orbit_y = a * np.sin(theta) * math.sin(math.radians(PLANETS[planet][2]))
# Draw orbit path (with transparency)
ax.plot(orbit_x, orbit_y, color='gray', alpha=0.3)
# Draw actual planet
size = 0.2 * max(0.3, min(1.0, np.log10(PLANET_SIZES[planet]) + 0.5))
planet_circle = Circle((plot_x, plot_y), size,
color=PLANET_COLORS[planet], zorder=5)
ax.add_artist(planet_circle)
# Label planets if requested
if show_names:
ax.annotate(planet, (plot_x, plot_y), xytext=(5, 5),
textcoords='offset points', fontsize=10, color='white')
# Add distance indicators (circular grid lines)
for distance in [5, 10, 15, 20, 25, 30]:
if distance <= max_au:
circle = plt.Circle((0, 0), distance, color="gray", fill=False, alpha=0.2, linestyle="--")
ax.add_artist(circle)
# Label the distance grid
ax.annotate(f"{distance} AU", (distance, 0), color="gray", alpha=0.7,
fontsize=8, ha='left', va='bottom')
# Set title and labels
date_str = date.strftime("%Y-%m-%d")
ax.set_title(f"Solar System on {date_str} ({view}-view)", fontsize=16, color='white')
ax.set_xlabel("x (AU)" if view == "top" else "x (AU)", color='white')
ax.set_ylabel("y (AU)" if view == "top" else "z (AU)", color='white')
# Grid
ax.grid(True, alpha=0.3)
# Remove ticks for cleaner look
ax.tick_params(axis='both', which='both', length=0, colors='white')
# Background
ax.set_facecolor('black')
return fig
def explore_solar_system(date_str=None, days_from_now=0):
"""
Explore the solar system on a specific date or days from today.
Parameters:
- date_str: Date string in 'YYYY-MM-DD' format
- days_from_now: Number of days from today (used if date_str is None)
"""
if date_str:
try:
date = datetime.strptime(date_str, "%Y-%m-%d")
except ValueError:
print("Invalid date format. Please use YYYY-MM-DD.")
return
else:
date = datetime.now() + timedelta(days=days_from_now)
# Create top and side views
fig_top = plot_solar_system(date, view="top")
fig_side = plot_solar_system(date, view="side")
plt.tight_layout()
plt.show()
# Print planetary positions
print(f"Planetary positions on {date.strftime('%Y-%m-%d')}:")
print("-" * 50)
print(f"{'Planet':<10} {'X (AU)':<10} {'Y (AU)':<10} {'Z (AU)':<10}")
print("-" * 50)
for planet in PLANETS:
x, y, z = calculate_planet_position(planet, date)
print(f"{planet:<10} {x:<10.3f} {y:<10.3f} {z:<10.3f}")
def animate_solar_system(start_date, days=365, interval=20):
"""
Create an animation of the solar system over a period of time.
Parameters:
- start_date: datetime object for the starting date
- days: number of days to animate
- interval: time between frames in milliseconds
"""
# Set up the figure
fig, ax = plt.subplots(figsize=(10, 10))
max_au = 32 # Maximum distance to show
# Function to initialize the plot
def init():
ax.clear()
ax.set_xlim(-max_au, max_au)
ax.set_ylim(-max_au, max_au)
ax.set_aspect('equal')
ax.set_facecolor('black')
return []
# Function to update the plot for each frame
def update(frame):
ax.clear()
ax.set_xlim(-max_au, max_au)
ax.set_ylim(-max_au, max_au)
ax.set_aspect('equal')
ax.set_facecolor('black')
# Current date
current_date = start_date + timedelta(days=frame)
date_str = current_date.strftime("%Y-%m-%d")
ax.set_title(f"Solar System on {date_str}", fontsize=16, color='white')
# Draw the Sun
sun = Circle((0, 0), 0.25, color='yellow', zorder=10)
ax.add_artist(sun)
# Draw planets and their orbits
objects = []
for planet in PLANETS:
# Get position
x, y, z = calculate_planet_position(planet, current_date)
# Draw orbit (simplified circular orbit)
a = PLANETS[planet][0]
theta = np.linspace(0, 2 * np.pi, 100)
orbit_x = a * np.cos(theta)
orbit_y = a * np.sin(theta)
ax.plot(orbit_x, orbit_y, color='gray', alpha=0.3)
# Draw planet
size = 0.2 * max(0.3, min(1.0, np.log10(PLANET_SIZES[planet]) + 0.5))
planet_circle = Circle((x, y), size, color=PLANET_COLORS[planet], zorder=5)
ax.add_artist(planet_circle)
objects.append(planet_circle)
# Label planet
ax.annotate(planet, (x, y), xytext=(5, 5),
textcoords='offset points', fontsize=10, color='white')
ax.grid(True, alpha=0.3)
ax.set_xlabel("x (AU)", color='white')
ax.set_ylabel("y (AU)", color='white')
ax.tick_params(axis='both', colors='white')
return objects
# Create the animation
ani = animation.FuncAnimation(fig, update, frames=days,
init_func=init, blit=False, interval=interval)
plt.tight_layout()
return ani
def main():
"""Main function to run the solar system visualisation program."""
print("Solar System visualisation")
print("=" * 30)
while True:
print("\nWhat would you like to do?")
print("1. View solar system on a specific date")
print("2. View solar system today")
print("3. Create animation")
print("4. Exit")
choice = input("\nEnter your choice (1-4): ")
if choice == '1':
date_str = input("Enter date (YYYY-MM-DD): ")
explore_solar_system(date_str=date_str)
elif choice == '2':
explore_solar_system()
elif choice == '3':
start_date_str = input("Enter start date (YYYY-MM-DD, default=today): ")
if start_date_str:
try:
start_date = datetime.strptime(start_date_str, "%Y-%m-%d")
except ValueError:
print("Invalid date format. Using today's date.")
start_date = datetime.now()
else:
start_date = datetime.now()
days_str = input("Enter number of days to animate (default=365): ")
days = int(days_str) if days_str else 365
print("Creating animation (this may take a while)...")
ani = animate_solar_system(start_date, days=days)
save = input("Save animation to file? (y/n, default=n): ")
if save.lower() == 'y':
filename = input("Enter filename (default=solar_system.mp4): ")
if not filename:
filename = "solar_system.mp4"
ani.save(filename, writer='ffmpeg', fps=30, dpi=100)
print(f"Animation saved to {filename}")
plt.show()
elif choice == '4':
break
else:
print("Invalid choice. Please try again.")
if __name__ == "__main__":
main()