HW3: The N-Body Problem (20 Points)

Overview / Logistics

The purpose of this assignment is to get you practice with the animation loop (i.e. change in time, update physics, update graphics) and nested loops in a fun, visual application. This assignment was heavily inspired by an assignment in Princeton COS 126: Intro To Computer Science.

What to submit: When you are finished, you should submit NBody.py to Canvas

The Problem

The video above shows how to create a simulation of the moon orbiting the earth using vpython (click here to see the accompanying code). In this assignment, you will extend this numerical integration scheme to approximate a simplified version of the n-body problem, in which n different point masses all exert gravity on each other at once. This scheme will not be entirely numerically accurate, especially since we are taking large steps in time during each iteration of the loop, but it will get the times of the orbits of each planet roughly correct, as you will see.

Net forces

One thing that makes this a bit trickier than just the earth and the moon is that for each body, all of the other bodies exert force on it at once. We refer to this as a net force. The acceleration due to all of these forces can simply be obtained at each timestep by summing the acceleration vectors for each individual body. The acceleration on body i due to body j, which we refer to as ai,j can be computed with the following formula

\[ a_{i,j} = \frac{G m_j}{||p_j - p_i||^3} (p_j - p_i)\]

where pi is the position of the ith rigid body, pj is the position of the jth rigid body, mj is the mass of the jth body, and G is the gravitational constant. To obtain the net acceleration at the ith body, we simply add all of the acceleration vectors for all bodies. For each planet i, this looks like

\[ a_i = a_{i,0} + a_{i, 1} + a_{i,2} + ... + a_{i,i-1} + a_{i, i+1} + ... + a_{i, N-1} \]

which can naturally be performed in a loop.

Code To Write

Click here to download the skeleton code for this assignment. You will be editing the file NBody.py. By default, this file loads the text file 4Planets.csv, which has information about the initial position, the initial velocity, and the mass of the Sun, Mercury, Venus, Earth, and Mars. (You are free to edit this file in Microsoft Excel or a similar program once you get this to work to see what happens when you change the masses and speeds). The code already loads the universe and sets up the VPython objects to display it. It also performs step 1 of the animation loop of computing elapsed time, as well as step 3 of updating the graphics. It will be your job to implement step 2 of the animation loop, which is updating the physics. There are three vpython vectors that you have to update during every iteration of the animation loop, in this order:

1. accelerations: This is a list of all of the acceleration vectors of the planets. Each acceleration starts off as (0, 0, 0) at each iteration of the loop. You will then have to update each vector each frame of the animation based on gravitation from the other planets.
2. velocities: This is a list of velocity vectors, one for each planet. You should update each one based on the accelerations you compute.
3. planets: This is list of vpython spheres, one for each planet. You should update the pos attribute of each of these based on its velocity, and this will update where the planets are drawn.

The pseudocode for the physics updates is as follows

Updating the acceleration

• For each rigid body i from 0 to n_planets-1
•            For each rigid body j from 0 to N-1, excluding i, compute the acceleration ai, j using the gravitation formula above, and add this to acceleration of planet i.

Updating the position/velocity

• For each rigid body i from 0 to N-1, update the ith row of V to be velocities[i] + dt*accelerations[i], and then update the ith planet's position to be planets[i].pos + dt*velocities[i].pos

If you have done this properly, you will see the following animation in your program

---