LAB 10


Planetary Motion

A scientist is tracking the position of a planet relative to a coordinate system whose origin is at the Sun and whose axes point at fixed stars. It is known that the orbit is an ellipse with the Sun at its focus. The major axis of the ellipse intersects with its circumference at two points known as the perihelion and aphelion:

It is also known that at the perihelion, the distance between the planet and the Sun is shortest and its orbital speed is largest. Whereas at at the aphelion, the distance between the planet and the Sun is largest and its orbital speed is smallest.

Beginning one day, to be called day #0, and about every week thereafter, the scientist records the day number (0, 1, 2, ...) and the x-y coordinates of the planet on a piece of paper. After several years, it was decided to automate the analysis of the recorded data and, hence, these pieces of paper were collected and entered, albeit not in order, into a text file named: orbit.txt (you can download this file by right-clicking its name and choosing to save it in the same directory in which you will write the program for this Lab). Each record in the file contains three real*8, space-delimitted numbers: the day number (not sorted), the x-coordinate, and the y-coordinate (measured in kilometers). For array declaration purposes, you can assume that the file has less than a thousand records.

We need to analyze this data to determine certain orbital attributes. To that end, we follow the structured approach and divide the problem into subproblems and then write the following main program:

	program orbit
	implicit none
	real*8 dayAR(1000), xAR(1000), yAR(1000)
	real*8 minMaxOrb(4), minMaxRad(4)
	real*8 min, max
	integer*4 lengthAR, maxSize /1000/, permute(1000)

	call loadFile(dayAR, xAR, yAR, lengthAR, maxSize)
	call sortMe(dayAR, xAR, yAR, permute, lengthAR)

	call radial(dayAR, xAR, yAR, lengthAR, minMaxRad)
	call orbital(dayAR, xAR, yAR, lengthAR, minMaxOrb, min, max)

	call showGeo(minMaxRad)
	call showGeo(minMaxOrb)
	call showPeriod(minMaxRad, min, max)

Here are descriptions (and hints) for each subprogram:


  1. Printed program
  2. A log showing the program's output (a single run).
  3. A brief explanation of how you computed the planet's period.
  4. This is optional: a quantitative explanation of the difference between the two methods of computing the semi-major axis length.