Deadlines

Initial deadline: Sunday Mar. 1, 23:59
Final deadline: Sunday Mar. 15, 23:59

Before you begin

Download the code skeleton for this lab, and check that you can compile and run it. (When run, you should see a window containing a black square with the heading "Newtonian Gravity Simulation", and a start/stop button.)
We recommend using IDEA for this lab. If you do, create a new project and put the files from the code skeleton in the src directory to get started.
Read this document in its entirety, and make sure you understand it. Ask a teacher if anything is unclear.

Learning outcomes

After completing this assignment, you will have:

created a real world physics simulation;
implemented a program with relatively less guidance than in the previous exercises;
practiced most of the material covered in the course so far;
practiced refactoring your code to improve encapsulation and abstraction.

Introduction

In this assignment, you will create a simulation of planets moving under Newtonian gravity. To make things simpler, the simulation only uses two dimensions.

Note: for this assignment you may want to create your own classes and/or write unit tests to help you structure your lab and ensure its correctness. If you do, please upload these files with your submission.

The Physics

Recall Newton’s Law of Gravity: between a point mass m₁ and a point mass m₂, there exists an attraction force directed along the line joining the point masses of size $LaTeX: \frac{Gm_1m_2}{r^2}$ where r is the distance between the point masses.

Units

Choosing the correct units makes building a simulation much easier. For this model, we will use the following units:

mass will be measured in solar masses, where 1 solar mass is the mass of the sun, approx. 2*10³⁰ kg;
distance will be measured in astronomical units (AU), where 1 AU is the average distance from the Sun to the Earth, approx. 1.5*10⁸ km;
time will be measured in years.

In these units, the constant G is approximately 4π² . (Can you see why? Hint: Assume the Earth’s orbit is perfectly circular. You do not have to answer this question as part of this lab, but it is a fun exercise.)

Part 1: The Planet Class

Write the class Planet. (Note: we are using the word ’planet’ for every body that we model. So, for this lab, we are calling the Sun a planet.) The class should have the following members:

instance variables giving the planet's position and velocity (in two dimensions);
an instance variable mass;
instance variables int radius and Color color
(these will be used only for displaying it on the screen, they do not affect the physics; Color is a class in the java.awt package);
a constructor that takes values for all the instance variables as parameters and sets them;
getters and setters for all the instance variables.

This is only our first version of the class - we will be modifying it later.

Part 2: The View

The next thing we are going to do is connect up the PlanetsView class, so that we can see the changes we make in future steps immediately.

Look at PlanetsModel. Add to this class the following:

Create two planet objects: sun and earth.
Put the Sun at coordinates (0,0) and stationary.
Put the Earth at coordinates (1, 0). It should initially be travelling vertically upwards with speed 2π (can you see why?).
The mass of the Sun is about 332946 times the mass of the Earth.
The update method is called 100 times a second when the program is run. Every time it is called, it should advance the model by dt years. For now, make it do just one thing when called: display the planets.
Replace the line // TODO Draw the planets with code that calls view.showPlanet() on both planets.

Run the program. You should see the Sun and the Earth stationary on the screen. Make sure this part is working before going any further!

Part 3: Euler's Method

We are going to use Euler’s method to simulate the motion of the planets first. Every time the update method is called, it should do the following:

calculate the acceleration of each planet, using Newton’s law of gravity and the position and mass of the other planet;
calculate the new position of each planet, using the formula: p(t + dt) = p(t) + v(t)dt;
calculate the new velocity of each planet, using the formula: v(t + dt) = v(t) + a(t)dt.

Write code that does this to replace the line // TODO Calculate the new positions of the planets. Run the program, and you should see the Earth travelling in a circle around the Sun. Again, get this part working before going any further.

Tip: you may have noticed that there's a lot of vector math creeping in here. This might be a good place to introduce a helper class or two of your own...

Part 4: The Three Body Problem

Add a third body to your model. You can make this anything you want, but for some suggestions:

Jupiter has a mass of about 0.001 solar masses, and is 5.2 AU from the Sun.
A typical asteroid has a mass of 10⁻⁸ solar masses, and is 3 AU from the Sun.
Halley’s comet has a mass of 10⁻¹⁶ solar masses, and is 0.58 AU from the Sun at its closest approach.

Try to get your code into a state where it would be possible to add a fourth body by adding/changing just 1–2 lines.

Try different initial velocities, and see what happens. (You may want to try to get realistic velocities, either by Googling, calculation or experimenting with your program.)

Part 5: Refactoring

This is good so far, but could be made better. So far, our Planet class is just a data class, and the PlanetsModel class is doing all the work, in danger of becoming a God class. It would be better to use the principle of Tell Don’t Ask: PlanetsModel should tell each Planet what to do, instead of asking for its position, velocity, etc.

So, add to the Planet class the following methods:

updateVelocity(double ax, double ay, double dt) - updates the planet’s velocity assuming that it travels at constant acceleration (ax, ay) over time dt;
updatePosition(double dt) - updates the planet’s position assuming that it travels with its current velocity over time dt.

Now, rewrite your update method in PlanetsModel using these methods. When you are done, you should be able to delete the methods getVX and
getVY from Planet.

Part 6: The Verlet Algorithm

The Euler method is not very accurate. Here is a more accurate method, called the Verlet algorithm or the second Taylor approximation. It finds each planet’s acceleration at the start and end of the time step, (i.e. times t and t + dt).

Let n be the number of planets. Let p_i(t) be the position of planet i at time
t; v_i(t) its velocity; and a_i(t) its acceleration (i = 1, . . . , n).

At the start of the main loop, we know p_i(t) and v_i(t) for all i.
Calculate the acceleration a_i(t) for all the planets, by applying Newton’s law of gravity using the positions p₁(t), . . . , p_n(t).
Estimate the position of the planets at the end of the time step using the
formula: $LaTeX: p_i\left(t+dt\right)=p_i\left(t\right)+v_i\left(t\right)dt+\frac{a_i\left(t\right)\left(dt\right)^2}{2}$
Calculate the acceleration a_i(t + dt) by applying Newton’s Law of Gravity to the positions p₁(t + dt), . . . , p_n(t + dt). (Note: you need to know all these positions before you can find a_i(t + dt).)
Estimate the velocity halfway through the time step by $LaTeX: v_i\left(t+\frac{dt}{2}\right)=v_i\left(t\right)+\frac{a_i\left(t\right)dt}{2}$
Estimate the velocity at the end of the time step by $LaTeX: v_i\left(t+dt\right)=v_i\left(t+\frac{dt}{2}\right)+\frac{a_i\left(t+dt\right)dt}{2}$

Now, play with your simulation! You can try simulating a binary star system, a star in orbit around a black hole, a probe traveling from Earth to Mars... See if you can observe this sort of chaotic behaviour, that shows why the three-body problem was so difficult to solve using algebra alone: https://vimeo.com/11993047.

Note: you may find that changing the type signatures of some methods will make this part easier. If you do, please document this in your submission. Of course, any such changes must not break the view or controller.

Part 7: Optional Gold-Plating

The simulation is more accurate if we slow it down (that is, decrease the time step) when it is changing rapidly. This will be when acceleration is high. Modify your code so that the time step decreases with the maximum of the accelerations of the planets.

Deliverables

All Java files you've created or modified. This includes:
- PlanetsModel.java
- Planet.java
- Any tests, help classes, etc. you created.
A README file briefly explaining your solution, as well as the purpose of any helper or test classes included in your submission. If you made changes to the types of any methods mentioned in this document (such as updateVelocity or updatePosition), please describe how and why.

For details about how to structure, package and submit your solution, see the general lab guidelines.