Physics 270, Problem 11: FORTRAN program for Kepler's motion

Physics 270 Problem 11: FORTRAN program to compute planetary motion (Kepler's laws)

1. Introduction

Write a FORTRAN program to compute the motion of the earth around the sun, assuming that the sun is located at the position (x,y)=(0,0). Use cartesian (x,y) coordinates.

The equations of motion are the following, in Euler's method:

r= sqrt(x*x + y*y)

accel_x = -GM*x/r*r*r
accel_y = -GM*y/r*r*r

vel_x(n+1) = vel_x(n) + accel_x(n)*dt
vel_y(n_1) = vel_y(n) + accel_y(n)*dt

x(n+1) = x(n) + vel_x(n)*dt
y(n+1) = y(n) + vel_y(n)*dt

2. Tasks

Use your program to calculate the trajectory for the earth under various initial conditions. Determine how to make the earth's motion circular, elliptical, and unbounded (doesn't return to the starting point).

Your program must save the trajectory as a function of time in a file. Your report must include a graphical plot of the trajectory for each of the three cases, along with the initial conditions. Try to express your results in terms of the parameters G and M.

In astronomical units, the product of G and M is a reasonable number. The unit of length in this system of units is called astronomical units, AU, and is based on the motion of the earth around the sun. The unit of length is the size of the semi-major axis of the elliptical orbit around the sun, a, which is defined in AU to be unity,

1 AU = 1.496 x 10^11 m = a The unit of time is taken to be one year, 3.15 x 10^7 s, the period of the earth's orbit around the sun. From Kepler's laws, the relationship between GM and the orbital motion of the earth is GM = (4pi^2*a^3)/T^2 = 4pi^2, in units of AU^3/yr^2 .