#!/usr/bin/python3 # =================================================================== # From: medium.com/analytics-vidhya/simulating-the-solar-system- # with-under-100-lines-of-python-code-5c53b3039fc6 # ---- # Note: The original code has some problems because some of the # syntax was deprecated and is no longer supported # (my changes are marked and shown in red) # ---- # Doc: astroquery.jplhorizons # astroquery.readthedocs.io/en/latest/jplhorizons/jplhorizons.html # ---- # Doc: astropy.time # docs.astropy.org/en/stable/time/index.html # =================================================================== import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation from astropy.time import Time from astroquery.jplhorizons import Horizons # ------------------------------------------------------------------- # ---- runtime # ------------------------------------------------------------------- # ---- simulating a solar system starting from this date sim_start_date = "2018-01-01" # ---- (int) simulation duration in days sim_duration = 2 * 365 # ---- mass of Earth relative to mass of the sun m_earth = 5.9722e24 / 1.98847e30 m_moon = 7.3477e22 / 1.98847e30 # ------------------------------------------------------------------- # ---- class for objects in the solar system # ---- define the objects: the Sun, Earth, Mercury, etc # ------------------------------------------------------------------- class Object: def __init__(self, name, rad, color, r, v): self.name = name self.r = np.array(r, dtype=np.float64) self.v = np.array(v, dtype=np.float64) self.xs = [] self.ys = [] self.plot = ax.scatter(r[0], r[1], color=color, s=rad**2, edgecolors=None, zorder=10) self.line, = ax.plot([], [], color=color, linewidth=1.4) # ------------------------------------------------------------------- # ---- class the solar system # ------------------------------------------------------------------- class SolarSystem: def __init__(self, thesun): self.thesun = thesun self.planets = [] self.time = None self.timestamp = ax.text(.03, .94, 'Date: ', color='w', transform=ax.transAxes, fontsize='x-large') # ---- add a planet to the solar system def add_planet(self, planet): self.planets.append(planet) # ---- evolve the trajectories def evolve(self): dt = 1.0 self.time += dt plots = [] lines = [] for p in self.planets: p.r += p.v * dt # acc is in units of AU/day^2 acc = -2.959e-4 * p.r / np.sum(p.r**2)**(3./2) p.v += acc * dt p.xs.append(p.r[0]) p.ys.append(p.r[1]) p.plot.set_offsets(p.r[:2]) p.line.set_xdata(p.xs) p.line.set_ydata(p.ys) plots.append(p.plot) lines.append(p.line) iso_time = Time(self.time, format='jd').iso # fixed code self.timestamp.set_text('Date: ' + iso_time) # fixed code return plots + lines + [self.timestamp] # ------------------------------------------------------------------- # ---- main # ------------------------------------------------------------------- # ---- amimate def animate(i): return ss.evolve() # ---- plot plt.style.use('dark_background') fig = plt.figure(figsize=[6, 6]) ax = plt.axes([0., 0., 1., 1.], xlim=(-1.8, 1.8), ylim=(-1.8, 1.8)) ax.set_aspect('equal') ax.axis('off') # ---- the solar system sun ss = SolarSystem(Object("Sun", 28, 'red', [0, 0, 0], [0, 0, 0])) ss.time = Time(sim_start_date).jd colors = ['gray', 'orange', 'blue', 'chocolate'] sizes = [0.38, 0.95, 1., 0.53] names = ['Mercury', 'Venus', 'Earth', 'Mars'] texty = [.47, .73, 1, 1.5] # ---- add the 1st, 2nd, 3rd, and 4th planets to the solar system for i, nasaid in enumerate([1, 2, 3, 4]): obj = Horizons(id=nasaid, location="@sun", epochs=ss.time, id_type=None).vectors() # fixed code: id_type ss.add_planet(Object(nasaid, 20 * sizes[i], colors[i], [np.double(obj[xi]) for xi in ['x', 'y', 'z']], [np.double(obj[vxi]) for vxi in ['vx', 'vy', 'vz']])) ax.text(0, - (texty[i] + 0.1), names[i], color=colors[i], zorder=1000, ha='center', fontsize='large') # ---- do the animation ani = animation.FuncAnimation(fig, animate, repeat=False, frames=sim_duration, blit=True, interval=20) plt.show() # ---- save animation to a file # ani.save('solar_system_6in_150dpi.mp4', fps=60, dpi=150)