Homework 8

Due: Tuesday, April 22 at 11:59 pm

Instructions

Please read the instructions carefully. An incorrectly formatted submission may not be counted.

There are several questions in this assignment. This assignment starts from skeleton code called nbody.py. Please include comments in your code detailing your thought process where appropriate. Put nbody.py in a folder called uni-hw8 (so for me this would be tkp2108-hw8). Then compress this folder into a .zip file. For most operating systems, you should be able to right click and have a “compress” option in the context menu. This should create a file called tkp2108-hw8.zip (with your uni). Submit this on courseworks.

N Body Problem

In this homework, we will continue the work started in HW7.

When we left off, we had all of the base calculations we needed to implement our simulation. Recall from class, when we built a 1-D simulation of the kinematic equations:

def run_sim(starting_pos, starting_v, delta_t):
    time = [0]
    position = [starting_pos]
    velocity = [starting_v]
    current_position = position[-1]
    current_velocity = velocity[-1]
    while current_position >= 0:
        new_position = pos_next(current_position, current_velocity, delta_t)
        new_velocity = vel_next(current_position, current_velocity, delta_t)

        position.append(round(new_position, 4))
        velocity.append(round(new_velocity, 4))

        current_position = new_position
        current_velocity = velocity[-1]

        time.append(round(time[-1] + delta_t, 4))
    return time, position

First, we started with initial positions and velocities. In each time step, we:

calculated the next position given the current velocity
calculated the next velocity given the current acceleration

For our n-body problem, we will do the same exact thing, just in 2 dimensions:

current positions and velocities are stored in our objects array
next x and y position is calculated as a function of current x and y velocity
next x and y velocity is calculated as a function of x and y acceleration, where:
a = F/m
F is the sum of forces in the x/y direction
m is the mass of the individual object

Problems

Complete the skeleton code below.

Code for animating is provided, as is a 3-part main function:

animate a 2-body system, which should oscillate
animate an n-body system with an initial rotation
(commented) animate the n-body system to a gif file instead of live

Important modifications

We’ve made gravity strong and our objects relatively large, so we need to make an important adjustment. In your calculation of the gravitational force, limit the distinance to 10 (our minimum object size), e.g. np.max([s, 10.0]) ** 2

Additionally, we need to give our objects an initial velocity in a rotating fashion. Modify your generateObjects to do so. As an example, assuming a 5xN array data:

def generateObjects(count, center=0, width=10, include_blackhole=False):

   <your code from hw7>

   if include_blackhole:
        # now calculate the center of mass
        com = centerOfMass(data)

        # now get the angles between the com and the data points
        angles = [angle(com, datum[:2]) for datum in data]

        # if each point is θ from com, then 90+θ (π/2+θ rad) should be good
        angles = angles + np.ones(count) * np.pi / 2

        # v = root(GM/r)
        velocity_magnitudes = np.sqrt(np.array([vec[2] / distance(com, vec[:2]) for vec in data]))

        velocities = np.array(
            [
                [np.cos(theta) * mag, np.sin(theta) * mag]
                for theta, mag in zip(angles, velocity_magnitudes)
            ]
        )

        # slice velocities to data
        data[:, 3:5] = velocities

        # as a last step, put a stationary black hole at the center
        data[0] = np.array([center, center, data[:, 2].max() * 1e4, 0, 0])
    return data

Notes

Due to an annoying bug in matplotlib/spyder, the animations may not render properly. To fix, do the following:

In Spyder settings, update the graphics backend to be “automatic”:

Demo

Full Skeleton Code

Full skeleton code from the solutions is available below.

NOTE: You are welcome to use your own code from HW7.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter


def animateObjectsFrame(objects, ax, lim=None):
    """Animate a matplotlib chart of our objects. `objects` argument is an Nx2 matrix of (x,y) coordinates"""
    # clear the frame
    ax.clear()

    # scatter the data, use max size for object size
    sizes = [min(100, x) for x in objects[:, 2]]
    ax.scatter(objects[:, 0], objects[:, 1], s=sizes, alpha=0.5)

    # plot the center of mass in red
    com = centerOfMass(objects)
    ax.scatter(com[0], com[1], s=25, c="red")

    if lim:
        ax.set_xlim(lim)
        ax.set_ylim(lim)

_animations = []

def animate(objects, time, delta_t, lim=None, save_gif=False):
    _animations.clear()
    frames = int(time / delta_t)
    fig, ax = plt.subplots(figsize=(14, 8))

    def animateStep(i):
        simulateStep(objects, delta_t)
        animateObjectsFrame(objects, ax, lim)

    ani = FuncAnimation(fig, animateStep, interval=100, frames=frames)
    _animations.append(ani)

    if save_gif:
        # Save as GIF
        writergif = PillowWriter(fps=30)
        ani.save("plot.gif", writer=writergif)
    else:
        # Show the plot
        plt.show()


def nextPosition(vec, delta_t):
    """Calculate the next position of `vec` given forces `forcevec`

    Args:
        vec: 1x5 vector of (x,y,mass,vx,vy)

    Returns: (new x position, new y position)
    """
    # TODO!
    return 0, 0


def nextVelocity(vec, forcevec, delta_t):
    """Calculate the next position of `vec` given forces `forcevec`

    Args:
        vec: 1x5 vector of (x,y,mass,vx,vy)
        forcevec: 1x2 vector of (Fx, Fy)

    Returns: (new x velocity, new y velocity)
    """
    # TODO!
    return 0, 0


def nextPositionAndVelocity(vec, forcevec, delta_t):
    """Calculate the next velocity and position of vector based on the forces acting on it"""
    # TODO!
    # NOTE: update `vec` in-place, don't return anything

def simulateStep(objects, delta_t):
    """Given the objects, which have the form [x, y, mass, xvel, yvel], calculate the sum of forces
    acting on each object and update the object with their next position and velocity"""
    # TODO!
    # NOTE: update `objects` in-place, don't return anything


def main():
    time = 1000
    delta_t = 10

    # testing with a binary system with 0 net initial velocity is a good bet
    data = generateObjects(2, width=10, center=0, include_blackhole=False)
    animate(data, time, delta_t, lim=(-100, 100))

    # now lets test with a realistic galaxy with a black hole
    # (note that i picked specific numbers here and played with the function's default arguments
    # in order to generate pretty picture. feel free to play around with your own values)
    data = generateObjects(100, width=1e5, include_blackhole=True)

    animate(data, time, delta_t, lim=(-1e3, 1e3))
    # To save a gif, commment the above and uncommment this
    # animate(data, time, delta_t, lim=(-1e3, 1e3), save_gif=True)


if __name__ == "__main__":
    main()