Assignment 8: Solar system

Due: 9:00am, Monday, March 8. Value: 30 pts. Submit to Sauron.

In this assignment we'll simulate the law of gravity to illustrate the motions of the sun and the planets of inner solar system: Mercury (cyan), Venus (green), Earth (blue), and Mars (red), as well as a comet based on Comet Encke (magenta), which has an orbit of 3.3 earth years.

The program is already started below; your job will be to complete the stepBody and updateBodyVelocity functions.

from graphics import *
import time
import math

steps_per_frame = 1  # increase this to make simulation go faster

def stepBody(body):
    pass  # replace this line with your implementation

def updateBodyVelocity(body, system):
    pass  # replace this line with your implementation

def createSystem(window, system):
    addBody(system, 1.9891e30,   0.000, 0.0000, 'yellow', 'Sun')
    addBody(system, 3.3022e23,  57.909, 0.4787, 'cyan',   'Mercury')
    addBody(system, 4.8685e24, 108.209, 0.3502, 'green',  'Venus')
    addBody(system, 5.9736e24, 149.598, 0.2978, 'blue',   'Earth')
    addBody(system, 6.4185e23, 227.939, 0.2408, 'red',    'Mars')
    addBody(system, 1.0000e20,  49.397, 0.7051, 'magenta','Encke') # Comet Encke

    for body in system:
        body[0].draw(window)

def addBody(system, mass, radius, vel, color, name):
    center = Point(250, 250 - radius)
    circle = Circle(center, 5)
    circle.setFill(color)
    body = [circle, vel, 0.0, mass, name]
    system.append(body)

window = GraphWin('Solar System', 500, 500)
system = []
createSystem(window, system)
while not window.isClosed():
    for step in range(steps_per_frame):
        for body in system:
            stepBody(body)
        for body in system:
            updateBodyVelocity(body, system)
    time.sleep(0.03)

Part 1: stepBody

Complete the stepBody function. This function takes a single parameter body, which is a list containing information about a body in the solar system:

body[0]:	A Circle as defined in the graphics module, which has been drawn into the GraphWin.
body[1]:	A number indicating the body's velocity in the x-direction.
body[2]:	A number indicating the body's velocity in the y-direction.
body[3]:	A number indicating the body's mass.
body[4]:	The name of the body. (You won't actually need to use this.)

The stepBody function's job is to tell the Circle to move according to the body's current velocity.

You should complete and test this function before moving to Part 2. When testing, you should find that all the bodies except the sun move slowly to the right (east), with the bodies closer to the sun moving more quickly. They don't go around the sun because you haven't yet added any gravitational effects.

Part 2: updateBodyVelocity

Now you should complete the updateBodyVelocity function, which updates a single body's velocity based on the gravitational pull exerted by the other bodies. The function takes two parameters: The first (body) contains information about a body whose velocity should be updated based on the gravitational pull by other bodies in the system, and the second (system) is a list of information about all bodies in the system (including body). Both body and each element in the system list are represented themselves as lists, as described in Part 1.

Note that this function should alter the velocity only of body, not of other bodies in the system. The program calls updateBodyVelocity for every body in the system, and each function call will update the velocity only for that one body.

(Yes, it will call updateBodyVelocity for the sun. Your program should not attempt to treat the sun specially: After all, the sun is just another body, which is pulled by the other bodies — it's just that its mass is so large that these pulls are negligible, and after half a year the planet has basically canceled out its pull by pulling the sun in the opposite direction anyway.)

So how do you compute the gravitational pull? Suppose we have two bodies, one of mass m₀ at (x₀, y₀), and another of mass m₁ at (x₁, y₁). We can determine the distance and angle between them as diagrammed below.

I recommend using the math module's atan2 function to compute θ. It takes two parameters, the first being the numerator of the above fraction, y₁ − y₀, and the second being the denominator, x₁ − x₀. It returns the angle in radians.

Gravity will lead this mass m₀ to be pulled toward m₁; the law of gravity says that it will be pulled with a force of:

(To simplify the discussion, I am calling this force, but technically this is the equation for the acceleration on the body created by the gravitational force.)

The gravitational constant G used in the above formula is 6.67 × 10⁻³⁰. (The actual gravitational constant is 6.67 × 10⁻¹¹ m³ / kg s²; but in our program each pixel represents 10⁹ meters, and each frame represents 10⁴ seconds; after converting the units to pixels and frames, we arrive at the exponent of −30.)

Of course, m₀ will already have some current velocity, moving v_x pixels in the x-direction with each time frame and v_y pixels in the y-direction. But the gravitational acceleration will lead v_x to increase by F cos θ and v_y to increase by F sin θ. The math module contains cos and sin functions that do this (each taking the parameter in radians, which conveniently is what atan2 returns).

Of course, there are several bodies in the solar system. The same body m₀ will be pulled by all other bodies; thus, to update the velocity of a body, you would go through all other bodies, compute the force pulling it toward that body, and update its velocity accordingly.