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?
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.
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.
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!
