CODEWITHSUNDEEP
pythonastronomyastroquery
atin
atin

Reputation: 893

Planet not following orbit in Solar System Simulation

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

Upvotes: 0

Views: 169

Answers (1)

KEIRSTEN C
KEIRSTEN C

Reputation: 96

I think it is because the position and velocity of Mercury are given in astronomical units and AU per day and the gravitational constant and the masses of the celestial bodies are given in SI units. It is probably messing up the calculation. Try converting all of them to SI units. Convert AU to meters and velocity from AU/day to m/s.

Upvotes: 2

Related Questions