# Computerised Solar System Simulation [closed]

In school I am doing a project where using python I simulate the solar system

The physics I know so far is newtons laws and Kepler’s laws and I have done a lot of research but there appears to be lots of different ways to model the solar system and often include very complex mathematics. I know it won’t be simple but I’m a bit lost at the moment and was wondering where I should go next and how I simulate this?

## closed as too broad by tpg2114♦, ZeroTheHero, sammy gerbil, Sebastian Riese, PulsarMay 10 '18 at 14:16

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

• NBabel might be a good start. There's a few Q&A here on different aspects of these simulations (I've answered a few, but another user, David Hammen, has some really good answers too). – Kyle Kanos May 8 '18 at 19:30
• Not much we can do without more info. I can post an answer with some suggestions about how I'd implement a simulation (with python code, even!) but the more accurate your simulation gets, the more complex it gets, necessarily. Mercury, for example, has significant orbital differences compared to what Newton predicts, due to the intense gravity. – Jakob Lovern May 8 '18 at 19:33
• I can take a crack at answering this question, but it won't be until later this evening (US central time zone). – David Hammen May 8 '18 at 20:14
• @ Jakob Lovern - Mercury has incredibly small orbital differences compared to what Newtonian mechanics predicts, about 43 arcseconds per century. An arcsecond is a small quantity, less than one one millionth of a revolution. 43 of them per century is incredibly small. – David Hammen May 8 '18 at 20:18
• Would Computational Science be a better home for this question? – Qmechanic May 9 '18 at 4:11

For modeling any $n$-body system, you are essentially modeling two systems of equations: \begin{align} \mathbf v&=\frac{\mathrm d\mathbf x}{\mathrm dt} \tag{1} \\ \mathbf a&=f\left(\mathbf x\right)/m\tag{2} \end{align} where $f(\cdot)$ is the gravitational force: $$\mathbf F_i=m_iG\sum_jm_j\frac{\mathbf r_j -\mathbf r_i}{\vert\mathbf r_j-\mathbf r_i\vert^3}\tag{3}$$ (the force on object $i$ is the sum of the forces of all the other objects $j$ in the system considered).

When it comes to these systems, the most simple updating scheme (Euler method) is inappropriate because it does not conserve energy, so one needs to rely on what are called symplectic integrators. Probably the easiest one to implement would be the velocity verlet method, which is also briefly described in the NBabel page I linked to in the comments: \begin{align} x_{i+1}&=x_i+v_i\cdot dt+\frac{1}{2}a_i\cdot dt^2 \\ v_{i+1}&=v_i+\frac{1}{2}\left(a_i+a_{i+1}\right)\cdot dt \end{align} with the acceleration being computed via (2). Here, $i$ refers to the $i$th object. I mention this in another post (probably a few others, but this one sprung to mind).

The outline of the algorithm is,

initialize_grid()
while t < t_end:
for each object i:
a[i] = compute_accelerations(r)
for each object i:
r[i] = update_positions(r[i], v, a, dt)
for each object i:
a[i] = compute_accelerations(r)
for each object i:
v[i] = update_velocities(v[i], a[i], dt)
t = t + dt
(print diagnostics here?)
end do
print_final_positions()


The compute_accelerations function comes from (2) using the gravitational force in (3) and the update_positions and update_velocities are the Velocity Verlet equations. Initializing the objects might be a difficult part, but you can probably set them all at $\theta=0$ and integrate from there.

For larger number of bodies $n$, this algorithm becomes terrible inefficient, as it requires $n^2$ operations each step for computing the force, which is pretty slow for larger values of $n$. Some advanced methods can improve this, but it's likely not needed for your needs; what is written here is just fine for the solar system.

This is just a general outline of the physics ideas.

The absolute simplest thing you can do is just maintain the $(x,y,z)$ coordinates of each planet, and update using Newton's laws. Basically $a=F/m$, then $\Delta v=a\Delta t$, and $\Delta x=v\Delta t$. Lather, rinse, repeat.

Historically, people used spherical coordinates and implemented Kepler's laws directly. This is less likely to produce significant rounding errors that accumulate over long periods, because you're explicitly using the fact that the orbits are periodic. What is kind of a pain about this method is that it's not trivial to implement the function $\theta(t)$ that describes the angular motion of the planet around the sun for a given time. There is a transcendental equation involved, which can, e.g., be solved by Newton's method if you know calculus.

All of the above describes using Kepler's laws or Newton's laws (which are equivalent under the two-body approximation). Kepler's laws are an excellent approximation for the planets, especially if you use current data and don't try to project very far into the past or future. The earth-moon system is different, however. It isn't well approximated as a two-body system. It acts more like a three-body system: earth-sun-moon. Because this was considered an important problem historically, there are specialized methods for approximating this system that go back hundreds of years.

• Not worthy of a downvote, but this ($\Delta v = a\Delta t, \Delta x = v\Delta t$) is symplectic Euler. While not quite astoundingly awful (the category into which basic Euler falls), it is still quite awful. The rest of this answer is very good -- except it's not that hard to solve Kepler's problem ($M = E - e\sin E$, solve for the eccentric anomaly $E$ given the mean anomaly $M$ and the eccentricity $e$). Newton's method works fine so long as the eccentricity is not high (IIRC, Newton's method won't go awry if $e$ is less than about 2/3), which is the case of the planets in the solar system. – David Hammen May 8 '18 at 23:18
• @DavidHammen: She doesn't say much about what level of education she's at. If she's a high school student, this needs to be very simple. – Ben Crowell May 9 '18 at 3:32
• Coding an orbit sim using Euler integration is a disappointing exercise, IME. You need to use a tiny time step and stick with small eccentricities or your orbits simply won't close up. I like synchronized Leapfrog integration: it's symplectic, easy to code, and it can be understood as using an acceleration that's essentially an average of the acceleration over the time step. – PM 2Ring May 12 '18 at 23:26
• @BenCrowell: yes I'm a high school student (17) sorry for not mentioning that it would explain why i am so confused – lucy May 18 '18 at 20:18
• @PM2Ring i have researched Leapfrog integration abit but sorry for asking what is a Sympletic Integrator? – lucy May 18 '18 at 20:19

Define a parent 'body' object. It'll contain the basic physics to work as a moving object. In pseudopython (2.7) it would look a lot like:

class body:
def __init__(self, mass, velocity, position):
self.mass = mass
self.velocity = velocity
self.position = position
def move_step(self, time_scale):
#Basically, steps the object along the simulation.
self.position += self.velocity * time_scale
#Not actually that simple, but time_scale is used to change
#how accurate the simulation is.


Note that I'm just gonna assume that velocity and position are defined as vectors, so that we can easily get their x, y, or z components, or get the magnitude and direction. Writing the vector object will be left up to the reader. Note that I didn't include any references to gravity or other objects. I've found that object-to-object interaction is annoying to deal with, so it would be better to simply have an object that handles the simulation.

Essentially, it loops through each body in the simulation and calculates the forces on the body. Since $\Sigma F = m*a$, you can get the $\Delta V$ of the object and apply it to that body's velocity. After every body's velocity is updated, the simulator loops through again and triggers their move_step method.

Sounds good, right? All you need is to define all the gravities acting on an object, then sum them and update.

Not so fast. You gotta factor in relativity if you want increased accuracy. In general it doesn't change much, but when you get to the scale of the Sun and Mercury, it becomes noticeable. I'd implement a method (or a lambda) that calculates Lorentz factors, and a conditional that chooses whether to use Newtonian or relativistic math based on some fancy factor.

### Collisions

Ok, so I can't do this one without making things simpler. Stuff where their sizes are massively different is easy. larger.mass += smaller.mass; del smaller

After that, there's so much complexity that I'm not really sure how to proceed. In essence, objects need to break apart, and they need to have shape. That's hard to do. I'd recommend stealing a lot of stuff from the bpy package.

Gotta go, sorry!

• Lorentz factors!? Seriously?? – Kyle Kanos May 8 '18 at 21:11
• Sorry! I got excited. It's not necessary, I admit. – Jakob Lovern May 8 '18 at 21:13
• No, not at all necessary. This also doesn't really seem to cover at all what OP wants or needs. And what's with the "einstein punching newton" image that isn't actually an image? – Kyle Kanos May 8 '18 at 22:56
• Honestly, I'm a lot less impressed with this answer in retrospect. I guess I was trying to develop the OP's system from scratch, but not actually give all the programming outright. Feature creep must've gotten to me. As for the fake image, I guess I figured either someone would recommend one or I'd draw and upload it myself. – Jakob Lovern May 8 '18 at 23:02
• In short, this answer is total crap and needs to be reworked or removed. – Jakob Lovern May 8 '18 at 23:02