**Lab Assignment 2: Gravity simulator** # Checkpoint 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 (`earthmoon.py` and `solarsystem.py`). 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 System.py.) 4. some test code for all this. You can start with `earthmoon.py`, 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/s^2. (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. # 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. * A **body** is the generic term for a star, moon, or planet. A body has a position, a velocity, and a mass. * A **system** is the collection of bodies that we want to simulate. 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 $v_x = 0$ and $v_y = 0$. The Earth has mass of about $5.9736 \times 10^{24}$ 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 \times 10^8$ 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 \times 10^8$ 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 $v_x = 0$ and $v_y = 1022$. The mass of the moon is about $7.3477 \times 10^{22}$ 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 \cdot m_{\textrm{moon}} / {r^2}$$ where $G = 6.67384 \times 10^{-11}$ is the universal gravitational constant, $m_{\textrm{moon}}$ 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 $d_x$ is the *x* distance between the Earth and the moon, found by subtracting $x_{\textrm{moon}}$ – $x_{\textrm{Earth}}$, and $d_y$ is similarly the *y* distance, we can use these values to compute the *x* and *y* components $a_x$ and $a_y$ of the acceleration: $$a_x = a \; d_x / r$$ $$a_y = a \; d_y / r$$ Total distance $r$ between the two bodies used in above given formula can be computed as square root of $d_x^2 + d_y^2$. 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 $a_x$ 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. 4. Update the velocity of the body using $a_x$, $a_y$, 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 $a_x$ by the timestep and add that to $v_x$. Similarly, multiply $a_y$ by the timestep and add that to $v_y$. 5. Update the position of the body using $v_x$, $v_y$, and the timestep. That will be the position at which it's drawn at the next timestep. 6. 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 cs1lib.py to the project, download the file [earthmoon.py](earthmoon.py), and add it to the project. Read it. Here's a screenshot taken in the midst of the animation:  Of course, earthmoon.py 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: ~~~ python __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 earthmoon.py, 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 earthmoon.py, 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: ~~~ python __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 earthmoon.py. (Of course, `__init__` is not called directly by code in earthmoon.py, 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 earthmoon.py. 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: ~~~ python (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 `earthmoon.py`, the line ~~~ python earth_moon.draw(WINDOW_WIDTH / 2, WINDOW_HEIGHT / 2, PIXELS_PER_METER) ~~~ 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 *i*th 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 *i*th 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, solar.py, 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](http://nssdc.gsfc.nasa.gov/planetary/factsheet/) on the Internet with the data. The table does not contain the mass of the Sun; use $1.98892 \times 10^{30}$ 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/s^2. 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. # 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. 5. (? 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 `solar.py`. 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. `earthmoon.py` (even if you didn't change it) 2. `body.py` 3. `system.py` 4. `solar.py` 5. Screenshots of simulations from `earthmoon.py` and `solar.py` 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.