Building a Digital Orrery

An exercise in Object Oriented Programming

Version 0.1

It is your goal in this exercise to construct a Digital Orrery. An orrery is a mechanical model of the Solar System. Here, we will generalize this to anything that is mechanically similar to the solar system: a collection of things bound gravitationally.

(image: wikimedia)

By J. S. Oishi (Bates College)



In [ ]:

    
import numpy as np
import matplotlib.pyplot as plt

%matplotlib notebook

Problem 1) Building a basic set of objects

Our first task is to map our problem onto a set of objects that we instantiate (that is, make instances of) in order to solve our problem.

Let's outline the scope of our problem.

A solar system exists in a Universe; here we can ignore the gravitational perturbation on the Solar System from the rest of the Universe. Our model will consist of a small number of bodies containing mass. It might also contain bodies without mass, so called "test particles".

The problem to be solved numerically is the gravitational N-body problem,

$$\ddot{\mathbf{r}}_i = -G\sum_{i \ne j} \frac{m_j \mathbf{r}_{ij}}{r_{ij}^3},$$

where $\mathbf{r}_{ij} \equiv \mathbf{r_i} - \mathbf{r_j}$. This task itself can be broken into two components:

the force calculation
the ODE integrator to advance $\mathbf{r}_i$ and $\dot{\mathbf{r}}_i$ forward in time

Problem 1a

In disucssion with a classmate, sketch out a set of classes that you will need to complete this project. Don't worry about things like numerical integrators yet.

Also, sketch out interfaces (start with the constructor), but don't worry about writing code right now.

Once you're done, find me and I'll give you the minimal list of objects.



In [ ]:

    
class MyFirstClass():
    pass

class MySecondClass():
    pass

Problem 1b

Wire them up! Now that you have the list, try them out. Python makes use of duck typing, you should too. That is, if your object has a mass m, a position r and a velocity rdot, it is a Body.



In [ ]:

    
r0 = np.array([0,0,0])
rdot0 = np.array([0,0,0])



In [ ]:

    
b = Body(1,r0, rdot0)

Problem 2

Now, we code the numerical algorithms. We're going to do the most simple things possible: a brute force ("direct N-Body" if you're feeling fancy) force calculation, and a leapfrog time integrator.

The leapfrog scheme is an explicit, second order scheme given by

$$r_{i+1} = r_{i} + v_{i} \Delta t + \frac{\Delta t^2}{2} a_{i}$$$$v_{i+1} = v_{i} + \frac{\Delta t}{2} (a_{i} + a_{i+1}),$$

where $\Delta t$ is the time step (which we'll just keep constant), and the subscript refers to the iteration number $i$.

Note that this scheme requires a force update in between calculating $r_{i+1}$ and $v_{i+1}$.

Problem 2a

Write a method that implements the force integrator. Test it on simple cases:

two equal 1 $M_\odot$ objects in your universe, 1 AU apart
a $1\ M_\odot$ object and a $1\ M_{\oplus}$ object, 1 AU apart



In [ ]:

Problem 2b Write the leapfrog integration as a method in the Universe class. Test it on one particle with no force (what should it do?)



In [ ]:

Problem 2c

Wire it all up! Try a 3-body calculation of the Earth-Sun-Moon system. Try the Earth-Jupiter-Sun system!

Challenge Problem

Construct a visualization method for the Universe class
Read about the Fast Multipole Method (FMM) here and implement one for the force calculation



In [ ]:

    
# good luck!