HW4: The N-Body Problem (20 Points)

Due Friday 3/6/2020

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, along with answers to the following as a comment on Canvas:

Did you work with a buddy on this assignment? If so, who?
Are you using up any grace points to buy lateness days? If so, how many?
Approximately how many hours it took you to finish this assignment (I will not judge you for this at all...I am simply using it to gauge if the assignments are too easy or hard)
Your overall impression of the assignment. Did you love it, hate it, or were you neutral? One word answers are fine, but if you have any suggestions for the future let me know.
Any other concerns that you have. For instance, if you have a bug that you were unable to solve but you made progress, write that here. The more you articulate the problem the more partial credit you will receive (fine to leave this blank)

The Problem

We showed in class how to simulate the orbit of the earth around the moon in 3D using Euler steps (Click here to see the 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 timesetp by summing the acceleration vectors for each individual body. The acceleration on body i due to body j, which we refer to as a_i,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 p_i is the position of the ith rigid body, p_j is the position of the jth rigid body, m_j 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 arrays that you have to update during every iteration of the animation loop, in this order:

A: This is an N x 3 array of acceleration vectors. There is one (x, y, z) vector per row for each of N rigid bodies. Each rows starts off as (0, 0, 0) at each iteration of the loop. It will be your job to fill in the rows of this, as explained below
V: This is an N x 3 array of velocity vectors, one for each planet. You should update each row of these based on the accelerations you compute.
P: This is an N x 3 array of position vectors, one for each planet. You should update each row of these based on the velocities you compute. These are what are drawn at the end of each loop.

The pseudocode for the physics updates is as follows

Updating the acceleration

For each rigid body i from 0 to N-1
For each rigid body j from 0 to N-1, excluding i, compute the acceleration a_{i, j} using the gravitation formula above, and add this to the ith row of A.

Updating the position/velocity

For each rigid body i from 0 to N-1, update the ith row of V to be V[i, :] + dt*A[i, :], and then update the ith row of P to be P[i, :] + dt*V[i, :]

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

Tips

Remember that adding or subtracting vectors with numpy is as simple as the + and - operation.
In numpy, the magnitude of a vector v can be computed with
This is actually just the pythagorean theorem; it's taking the square root of the sum of the squares of components. Of course, you will have to substitute v for whatever vector of which you want the magnitude.
To compute the net acceleration, you will need to use a nested for loop. You'll need a for loop on the outside and a for loop on the inside. Below is an example of a nested for loop that prints five rows of 10 asterixes
```
**********
**********
**********
**********
**********
                                                
```
Be sure that you tab things properly! Everything inside the first loop (including the inner loop), should be tabbed over exactly once relative to the surrounding code, and everything inside the inner loop should be tabbed over exactly twice