Lab Assignment 2: Gravity simulator


Read the whole lab description first. Draw some pictures. Make sure you understand how the Body and System classes will work together with the simulators ( and Once you've done that, return to this section and read about what you'll need to do for the checkpoint.

For the checkpoint, implement basic Body and System classes. You should have:

  1. complete init methods for both classes, with all instance variable that you will need. (The Body class will have x, y, v_x, v_y instance variables, at least, to represent the location and velocity of the Body. It will probably not need acceleration instance variables, since you'll later write a method to compute those quantities at every time step. The System class should have least a list of Body objects as an instance variable.)

  2. draw methods for both classes

  3. update methods for the location and velocity of Body objects. You won't need any fancy physics for this -- the x and y locations are updated by the appropriate velocities times the time step, and the velocities are updated by some accelerations that will be passed into the method, as well as the timestep (also passed in). You'll be testing that this works by making up acceleration values, as described in the next bullet point. (Later on, you'll replace these made-up numbers by values computed in

  4. some test code for all this. You can start with, but will want to modify it. (Later, you can go back to the original version, once you've tested and submitted your code for the checkpoint.) Don't be afraid to get a little non-physical with your simulator for testing purposes. For example, I'd set a small acceleration in the positive x and y direction for the moon, and I'd expect it to move in a straight line up and to the right of my screen, with increasing speed. That's not how gravity works, but this will test your update methods to see if they are reasonable. How do you choose these numbers? I looked up the acceleration of the moon towards the Earth and found that it is about .0028 m/s2. (A previous version of these notes gave a much larger value, because it did not use the timescale constant.)

This completes the checkpoint. Once you've completed the checkpoint, the main challenge is to compute correct accelerations using the positions and masses of the bodies, as we discussed in lecture, and as is discussed in the rest of the lab. Of course, you'll also need to write the bigger solar system simulation code as well.

For Lab2 checkpoint you need to submit, and files along with a screenshot of your output.

Problem description

The goal of this lab is to design and implement a simulator that animates motion of stars, planets, or moons. I'll start by describing the problem in English. First, some definitions.

We have to decide on some units to measure the various quantities. Let's use meters (m), kilograms (kg), and seconds (sec). We'll implement the simulator in only two dimensions, so all positions and velocities will have only x and y components.

Let's look at an example, the Earth-moon system. The Earth-moon system has two bodies. Can you guess what they are? To start with, we'll put the Earth at the position x = 0, y = 0, although we expect the Earth to move slightly over time as the moon's gravitational force acts on the Earth. We'll also start the Earth standing still, so that its initial velocity components are vx = 0 and vy = 0. The Earth has mass of about 5.9736 × 1024 kg; I looked it up on the Internet.

Where is the moon? I found that the average distance between the center of the Earth and the moon is about 3.844 × 108 meters. Because the moon's distance from the Earth changes throughout its elliptical orbit, that's just a rough estimate, but it's good enough for our purposes. So let's start the moon at the position x = 3.844 × 108 and y = 0. We could start the moon at some other position, but placing the moon initially on the x axis will simplify the setup.

I found that the speed of the moon is about 1022 meters per second. Remember that velocity has a direction and a magnitude, with speed being the magnitude. The moon's velocity should be in a direction perpendicular to the line from the center of the Earth to the moon, since the orbit of the moon is fairly close to circular:

So I initialized the velocity of the moon to vx = 0 and vy = 1022. The mass of the moon is about 7.3477 × 1022 kg.

Mathematics and physics

Here is a review of some of the mathematics and physics involved in computing the motion of the Earth-moon system. We're going to repeatedly update the positions and velocities of the bodies in the system.

Here's how it will work:

  1. Draw each body at its current position.

  2. Compute the timestep that we would like to update the simulation by for each frame of the animation. Let's say that the frame rate for our simulator is 30 frames per second. But if we moved the simulation at the same rate as the world we live in, it would take a month for the moon to revolve around the Earth. That's well beyond the limit of our patience! So we need some factor by which to speed up the simulation. Let's call this factor the time scale. A day is 86,400 seconds long (60 seconds/minute × 60 minutes/hour × 24 hours/day). Suppose that the simulator simulates 100,000 seconds of movement for every second of real time. Then we could watch an entire day of simulated time in under one second. The moon takes approximately 27.3 days to make a full revolution around the Earth, and so we can see one revolution in about 23.6 seconds. The timestep to update the simulation in each frame of the animation would be 100,000 / 30, or about 3333.33 seconds. This will be the timestep for our simulation—the number used to update the velocity (using the acceleration) and to update the position (using the velocity).

  3. Compute the acceleration of each body. To compute the acceleration of the Earth caused by the moon, we can use the formula

a = G * mmoon/r2

where G = 6.67384 × 10–11 is the universal gravitational constant, mmoon is the mass of the moon, and r is the distance from the Earth to the moon.

We're not quite done. a is the total amount of the acceleration, but we need to know how much of the acceleration is in the x direction, and how much is in the y direction. If dx is the x distance between the Earth and the moon, found by subtracting xmoonxEarth, and dy is similarly the y distance, we can use these values to compute the x and y components ax and ay of the acceleration:

ax = adx/r ,
ay = ady/r

Total distance r between the two bodies used in above given formula can be computed as square root of dx2 + dy2.

The situation is slightly more complicated if there are more than two bodies in the system. To compute the total x acceleration of a body, we add up the ax contributions of all the other bodies. Let's take as an example the planets in the solar system. The total x acceleration of the Earth would be the x acceleration of the Earth due to the Sun, plus the x acceleration of the Earth due to Mercury, plus the x acceleration of the Earth due to Venus, etc.

  1. Update the velocity of the body using ax, ay, and the time step. If the timestep is small enough (and 3333.33 seconds is small enough for our purposes), we can pretend the accelerations are constant throughout the timestep, and multiply ax by the timestep and add that to vx. Similarly, multiply ay by the timestep and add that to vy.

  2. Update the position of the body using vx, vy, and the timestep. That will be the position at which it's drawn at the next timestep.

  3. Repeat starting at step 2.

Design of the Python program

The simulation would be easy to write if we had custom data types to define a body and a system. You could then use methods of the Body and System classes that update positions and velocities using the equations described above and that draw the System and Body. Your first job is to write the System and Body classes.

To get you started, I've written the code that contains the initial positions and velocities of the Earth and moon, creates Body objects for the Earth and moon, and adds these Body objects to a new System. We then perform the animation.

I recommend that you create a new project for this lab assignment. Add to the project, download the file, and add it to the project. Read it.

Here's a screenshot taken in the midst of the animation:

Of course, won't work until you've written the Body and System classes.

The Body class

The Body class represents an individual body, such as the Earth, the moon, the Sun, or any other planet. The version I wrote had four methods:

__init__(self, mass, x, y, vx, vy, pixel_radius, r, g, b)
update_position(self, timestep)
update_velocity(self, ax, ay, timestep)
draw(self, cx, cy, pixels_per_meter)

The first five parameters (other than self) to the __init__ method (mass, x, y, vx, and vy) describe the physical properties of the body at the start of the simulation. The next parameter (pixel_radius) is the radius of the circle used to draw the body, in pixels. The last three parameters (r, g, and b) specify the color in which to draw the body. Although pixel_radius is in pixels, it is used only for drawing the planet and will not be used in the physics calculations in any way. I picked numbers that looked good for my program. For example, for a solar system simulation, I chose the radius of the Earth to be 6 pixels, and of Venus to be 5 pixels.

I will let you figure out the meanings of the parameters to update_position and update_velocity on your own.

The draw method draws the body on the screen, calling draw_circle to actually do the drawing. The parameters cx and cy represent the location, in pixels, of the center of the simulation. A body at the location x = 0, y = 0 should be drawn at the center of the window. Notice that in, I centered the Earth at x = 0, y = 0, but in order to get the Earth to be drawn in the center of the window, I supplied the values WINDOW_WIDTH / 2 and WINDOW_HEIGHT / 2 as cx and cy. The parameter pixels_per_meter gives the scale of the simulation, and it's used to convert the position of the body (stored in meters) into pixel coordinates in the window. In, I wanted 10,000,000 meters to scale to 3 pixels, and so the value for this parameter was 3 / 1e7.

You may use the method headers I wrote and write just the bodies of the methods.

The System class

The System class I wrote had just one instance variable, which was a reference to a list. Each item in the list is a reference to a Body object.

Code that uses the System class interacts with it through the following methods:

__init__(self, body_list)
update(self, timestep)
draw(self, cx, cy, pixels_per_meter)

The System class may have more methods than these, but it has to have at least these methods. You'll notice that they're called by code in (Of course, __init__ is not called directly by code in, but it is called by the System constructor.)

The draw method of System just calls the draw method on each body in the body list. Like the draw method, you can see a call of the update method in The update method computes the accelerations on each body and then calls methods in Body to update the velocity and position of each body.

I wrote another method, compute_acceleration, in System. This method is called only by update in System. The compute_acceleration method computes the x and y components of the acceleration of the body at index n in the list. (The method takes the index n as a parameter.) Remember that to compute the acceleration, you will have to loop over all other bodies to add up their contributions to the acceleration of body n.

A couple of notes about how I implemented compute_acceleration. First, when looping over the bodies, I had to make sure that I didn't compute the acceleration of n on itself. So the body of the loop has a check to skip over computing the acceleration of body n on itself. Second, compute_acceleration needs to return not one, but two values: the acceleration in x and the acceleration in y. This requirement presents a little problem: a function can return only one thing, and since a method is a function, a method can return only one thing. So how can I make compute_acceleration return two values?

It turns out that Python has a structure called a tuple, and you can form one just by putting some values in parentheses and separating them by commas. For example, here's a tuple containing the first four prime numbers: (2, 3, 5, 7). A function (or method) is allowed to return a tuple, because a tuple is one "thing." So you could have compute_acceleration return a tuple with two values, the accelerations in x and y.

But how do you save these returned values? You can save them into a tuple, and just pull out pieces of the tuple as necessary. For example, you could call compute_acceleration within update as follows:

(ax, ay) = self.compute_acceleration(n)

And then later on, you can use the variables ax and ay as usual. This tuple feature of Python is very nice for making functions (or methods) that return multiple values. I've often wished that other programming languages provided it.

Programming note: The draw methods

Did you notice that I've mentioned two methods, both named draw? One is in the Body class, and it draws a single body. The other is in the System class, and it draws the entire system by drawing each body in the system. When you make a method call and the method name is in more than one class, how do you know which method is being called? Look at the object reference to the left of the dot in the method call. What kind of object does it reference? Whatever kind of object it references gives the class that the method is in.

So, in, the line


must call the draw method of the System class, because earth_moon is a reference to a System object.

Programming note: Stringing together references

In some of the methods in my System class, I needed to say something like the following:

  1. Go to the list of bodies that the instance variable references.
  2. Find the ith reference to a body in that list.
  3. Get the x-coordinate of that body.

How do I write that in Python? Let's tease it out. Suppose that my instance variable, which references a list of references to Body objects, is body_list. Then:

  1. To refer to that list, I would write self.body_list.
  2. To refer to the ith reference in a body to that list, I would write self._body_list[i].
  3. To get the x-coordinate of that body, I would write self.body_list[i].x.

Notice that I've strung together three references. self references the System object I'm working on. self.body_list references a Python list. self.body_list[i] is a reference to a Body object, so that I can get x of this Body object by writing self.body_list[i].x.

Part two: The solar system

Once you have completely written and tested your program with the Earth-moon system, write another Python program,, that simulates the motion of the Sun and first four planets of the solar system. This should be the last step you take, and it should be similar to the Earth-moon system, except that you will need more bodies and with different values for masses, initial positions, and initial velocities.

I found a table on the Internet with the data. The table does not contain the mass of the Sun; use 1.98892 × 1030 kg.

Here's a screenshot taken in the midst of my solar system simulation:

How to get started

There's a lot to do; it's important to tackle things piece by piece. I would recommend thinking through the design first, seeing how everything will work, and then writing small pieces of the program that you can test. There's no point in writing the whole enchilada and then when it doesn't work, you don't know where to start looking for the bugs. Write a small piece, and test it extensively so that you're satisfied that it works correctly. Then write another small piece, possibly interacting with the first small piece, test this new small piece extensively, and so on. If you operate in this way, when something doesn't work, you'll know that the error is most likely in the newest code you added.

For example, let's say you just want to see the Earth and moon drawn in the window at their initial positions, but not moving. You could write a version of the System class for which the body of the update method is just the Python command pass, meaning to do nothing. Then you wouldn't need compute_acceleration yet, either, and the update methods of the Body class could also just be pass. Eventually you have to write those methods, but not yet.

The first things you would write would then be the __init__ methods for each class, and the draw methods for each class.

Once you have the drawing working, you could work more on the Body class, and write the two methods that update the position and velocity. Test them. For example, at first you could choose accelerations of 0, and make sure that the position update is working.

Once it's time to write the update method of System, take it one step at a time. Make sure that your loops loop over the correct things by adding print statements. Make sure you are computing distances and accelerations accurately. For example, an Internet search tells me that the magnitude of the acceleration of the moon toward the Earth is 0.0027 m/sec2. Is this the value you get upon the first time through the simulation? Print out the value and temporarily add a call to cs1_quit to your program to make sure the value you get matches this before moving on.


You will be graded on the following criteria:



Correctness counts for only 22 out of 40 points. If your program works perfectly but is stylistically a mess, your grade will be low.

Extra credit ideas

Any extra credit work should be written in separate Python files. Don't add extra credit work into your main solution code. These are just a few ideas; the first three might not be too difficult. You can make up your own ideas for extra credit.

  1. (10 points) Improve the look of the program by using images to draw the planets and Sun, calling the appropriate functions from cs1lib. NASA has a website with public-domain photographs of planets and the Sun. The filename of the image should probably be a parameter to the __init__ method of a body, instead of the radius and color parameters. Remember that the image file needs to be in the same project as your Python code for load_image to work.

  2. (5 points.) Add a new system, in a new Python file. One possibility is the Sirius AB binary star system. You can find some data online about the system, but be careful—the orbit is highly elliptical, so you'll need to make sure that the distance you initially choose between the stars matches the velocity. (If you can't find the velocity data, choose a distance, and then choose different velocities until the orbit has the right properties.)

  3. (10 points) Higher-precision simulation. The accuracy of the simulator depends on the timestep, which in turn depends on the frame rate. A smaller timestep means that the simulation will be more accurate. Why? In reality, the velocities and accelerations are not constant during the timestep. You'll get less error if the timesteps are smaller. If you want to see actual error with the timesteps I've chosen, try running your Earth-moon simulation for a few minutes. (Real minutes, not simulated minutes.) If your code is right, you'll notice that the Earth drifts toward the edge of the window.

    One way to improve the precision is to write your update method to subdivide the timestep it is given into many smaller timesteps, and add up the results of the smaller timesteps. For example, if the timestep given to update is 3333 seconds, your update method could compute 100 updates with timesteps of 33.33 seconds, accumulate the results, and return that value.

  4. (10 points) The midpoint method. The accuracy of the simulator also is limited by the x = x + t * vx and vx = vx + ax * t type of timestep. Because the acceleration is computed only at the beginning of the timestep, the error in the velocity computation causess the planets to drift slowly away from the Sun. You can improve the update calculations using the midpoint method. Compute the acceleration at the beginning of the timestep. Then advance the simulation by one timestep (using the current velocity and timestep to move the body), and compute the acceleration again. Take the average of the accelerations that you computed at the beginning and end of the timestep. Go back to the beginning of the timestep, and use that averaged acceleration to actually do the update and get the new velocity. You can use a similar procedure to update the position using the velocity.

  1. (? points) Write an educational game that allows planets to be added using the mouse. Perhaps the planets would grow more massive the longer the mouse button was held down? Or add a ship with thrusters. Be creative!

What to turn in

If you discussed the assignment with anyone, please note the name of that person, and the nature of that discussion, in the comments of Do not forget to put your own name and date in the comments in a header, indicating that you wrote the submitted code.

Turn in the following:

  1. (even if you didn't change it)
  5. Screenshots of simulations from and
  6. If you do any extra credit, submit it in separate files, clearly labeled.

Honor Code

The consequences of violating the Honor Code can be severe. Please always keep in mind the word and spirit of the Honor Code.