I'm trying to simulate a solar system, but the planets don't seem to follow the orbit, instead, the distance between a planet and the sun increases. I can't figure out what's wrong with this code. Here is an MRE.
import numpy as np
import matplotlib.pyplot as plt
GRAVITATIONAL_CONSTANT = 6.674e-11
EARTH_MASS = 5.972 * (10**24)
SUN_MASS = 332954.355179 * EARTH_MASS
MERCURY_MASS = 0.06 * EARTH_MASS
AU2M = 1.495979e11
AUD2MS = AU2M / 86400
class SolarSystemBody:
def __init__(self, name, mass, position, velocity):
self.name = name
self.mass = mass
self.position = np.array(position, dtype=float) * AU2M
self.velocity = np.array(velocity, dtype=float) * AUD2MS
self.acceleration = np.zeros(3, dtype=float)
def update_position(self, dt):
self.position += self.velocity * dt + 0.5 * self.acceleration * dt**2
def update_velocity(self, dt):
self.velocity += self.acceleration * dt
def gravitational_force(self, sun):
r = sun.position - self.position
distance = np.linalg.norm(r)
direction = r / distance
force_magnitude = GRAVITATIONAL_CONSTANT * self.mass * sun.mass / distance**2
return force_magnitude * direction
def calculate_acceleration(self, sun):
force = self.gravitational_force(sun)
self.acceleration = force / self.mass
mercury_position = [0.1269730114807624, 0.281031132701101, 0.01131924496141172]
mercury_velocity = [-0.03126433724097832, 0.01267637703164289, 0.00390363008183905]
sun = SolarSystemBody("Sun", SUN_MASS, [0, 0, 0], [0, 0, 0])
mercury = SolarSystemBody("Mercury", MERCURY_MASS, mercury_position, mercury_velocity)
dt = 3600 * 24
total_time = 365 * dt
pos = np.zeros((365, 3), dtype=float)
i = 0
for t in np.arange(0, total_time, dt):
print(np.linalg.norm(sun.position - mercury.position))
pos[i, :] = mercury.position
mercury.calculate_acceleration(sun)
mercury.update_velocity(dt)
mercury.update_position(dt)
i += 1
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot(pos[:, 0], pos[:, 1], pos[:, 2])
plt.show()
Getting the initial data using astroquery
's HorizonsClass
from astropy.time import Time
from astroquery.jplhorizons import Horizons
import numpy as np
sim_start_date = "2023-01-01"
data = []
planet_id = 199 # Mercury
obj = Horizons(id=planet_id, location="@sun", epochs=Time(sim_start_date).jd, id_type='id').vectors()
name = obj["targetname"].data[0].split('(')[0].strip()
r = [np.double(obj[xi]) for xi in ['x', 'y', 'z']]
v = [np.double(obj[vxi]) for vxi in ['vx', 'vy', 'vz']]
Update:
Apparently mercury was slingshotting using the sun thats why the sudden increase in it's distance to sun.The issue seems to be that dt
was too big so acceleration wasn't getting updated quickly. Decreasing dt fixes the issue.
dt = 360
total_time = 365 * 240 * dt
source https://stackoverflow.com/questions/76133729/planet-not-following-orbit-in-solar-system-simulation
Comments
Post a Comment