Project 8: The orbit of the Earth

Wednesday, Apr 30, at 5:00 PM: Neatly written equations showing EXACTLY how position and velocity will be calculated in the x- and y-directions
Monday, May 5, at 9:00 AM: The finished MATLAB code and analysis

Your job in this project: calculate the motion of the Earth around the Sun. It should move in an elliptical orbit ... but can you make it do so?

The situation

To a very good approximation, we can consider the Earth and Sun to be alone in the universe. The two objects attract each other via the force of gravity. Again to a good approximation, we may consider the Sun to be so much more massive than the Earth that it doesn't move appreciably. Thus, the only motion we need to follow is that of the Earth. We can further simplify things by considering only motions in the two-dimensional plane of the Earth's orbit.

You can use Euler's method for this project. The basic structure of the program will look something like this:

use the current position of the Earth to calculate the gravitational force on it
break the gravitational force into x and y components
use the force components to update the velocity of the Earth in each direction
use the velocity components to update the position of the Earth in each direction
update the time
go to Top of Loop

The tricky part will be figuring out how to break the gravitational force into its x- and y-components. Be sure to write down these equations explicitly and clearly in the pseudocode you submit on Wednesday.

You'll need some initial conditions. Use the following:

initial time is t=0 seconds, on midnight of Jan 1, 2000
the Sun is now and always fixed at the origin
at this time, the Earth is at position (x = 1.5e11 m, y = 0 m)
at this time, the Earth has velocity (vx = 0 m/s, vy = +29,747 m/s)

Adopt the following values:

mass of Sun is 1.99e30 kg
mass of Earth is 5.98e24 kg

The assignment

You should write a MATLAB routine which calculates the position and velocity of the Earth as a function of time. It should look like this:

  function orbit(timestep, duration, output_file)

where

     timestep                is the size of the timestep to use in 
	                           calculations, in days

     duration                is the length of time over which to simulate
                                   the orbit, in days

     output_file             is the name of a file in which you will
                                   write values of time, position, velocity

As you go through a loop, calculating position and velocity of the Earth, you should write the values into a text file, for later use.

The output file should contain data like this:

  0.0    1.5000000e+11 0.0000000e+00    0.0000000e+00 2.9747000e+04 
  1.0    1.5000000e+11 2.5701408e+09   -5.0969472e+02 2.9747000e+04 
  2.0    1.4995596e+11 5.1402816e+09   -1.0191651e+03 2.9738271e+04

where the first column is time (in days) since the start of the simulation, the second column is x-position (in meters), the third column is y-position (in meters), the fourth column is x-velocity (in meters per second), and the fifth column is y-velocity (in meters per second).

You must do the following:

Write a program which calculates the orbit of the Earth, as described above.
Run your program using a timestep of 1 day, and a duration of 20 days. Look at the results. Plot the position of the Earth, using a graph which runs from x = -1.8e11 to x = 1.8e11 meters, and y = -1.8e11 to y = 1.8e11 meters. Does everything look okay?
Calculate the total mechanical energy of the Earth-Sun system, by adding together the kinetic energy of the Earth and the gravitational potential energy of the Earth (relative to the Sun). Compare the total energy of the system at the start and end of your simulation. Write them down. They should be identical. Are they? If not, by what percentage do they differ?
Run your program using a timestep of 1 day, and a duration of 365 days. Again, examine and plot the results. Does everything look okay?
Again calculate the total mechanical energy of the Earth at the start and end of the simulation. Write them down. Are the values the same? If not, by what percentage do they differ?
Run your program with a timestep of 0.1 days. Calculate the change in total mechanical energy from start to 20 days, and from start to 365 days. How do the changes compare to your earlier values? Discuss.

Bells and Whistles

Try to run your simulation for many years. How many years can you simulate before the orbit goes bad? You might define "going bad" by a big change in total energy, or by a big change in shape or size of the orbit. Does it help to change the timestep?
Add the Earth's moon to the simulation, and make a plot showing the orbits of the Earth and Moon as they go around the Sun during the course of one year.
Show that your positions of the Earth define an ellipse.
Write a function which calculates the angular momentum of the Earth around the Sun at any given moment. Use it to track the angular momentum of the Earth in its orbit. Just like energy, angular momentum should be conserved. Is it?
Calculate the motion of the Sun due to the gravitational force of the Earth.

This page maintained by Michael Richmond. Last modified May 5, 2003.