                                Program G3BODY

              A Computer Graphics Display of the Three-Body Problem

                     Copyright (c) 1989-1994 by David Eagle


Program G3BODY is a computer program for the IBM (tm) and true compatible 
computers which can graphically display orbital motion in the circular 
restricted three-body problem (CRTBP). The software supports the VGA graphics
mode of the IBM-PC and true compatible computers.

This program first appeared in the Recreational Computing column of the Winter 
1989-90 (Volume 2, Number 1) issue of "Celestial Computing".

This algorithm and program examples are based on the methods described in 
"Periodic Orbits in the Restricted Three-Body Problem with Earth-Moon Masses", 
by R. A. Broucke, JPL TR 32-1168, 1968. Additional information can also be 
found in the JPL Technical Report, "Periodic Orbits in the Elliptic Restricted 
Three-Body Problem", by R. A. Broucke, JPL TR 32-1360, 1969. This report 
contains a listing of a FORTRAN computer program which also calculates 
periodic orbits in the "elliptic" restricted three-body problem.

Program G3BODY simulates the motion of a spacecraft in the Earth-Moon system. 
In the discussion which follows, subscript 1 is the Earth and subscript 2 
represents the Moon.

The motion of the spacecraft is displayed in a coordinate system which is 
rotating about the center-of-mass or barycenter of the Earth-Moon system. The 
orbital motion is confined entirely to the x-y plane.

For convenience the problem is formulated in "canonical" units. The unit of 
length is taken to be the (assumed) constant distance between the Earth and 
Moon, and the unit of time is chosen such that the Earth and Moon have an 
angular velocity w about their barycenter equal to one. With this formulation, 
Kepler's third law is as follows:

	 w^2*abs(m1*m2)^3 = g(m1 + m2) = 1            

and the value for the universal gravitational constant g is 1.

The x and y coordinates of the Earth and Moon in this coordinate system are 
as follows:

	 x1 = -mu,     y1 = 0,     x2 = 1 - mu,     y2 = 0.

In his JPL technical reports, Professor R. Broucke calls these "synodical" 
coordinates.

Program G3BODY numerically integrates the first order form of the vector 
differential equations of motion using a fourth-order, variable-step-size 
Runge-Kutta-Fehlberg method.

In his prize memoir of 1772, Joseph-Louis Lagrange proved that there are five 
equilibrium points in the restricted three-body problem. If we place a 
satellite at one of these points with zero initial velocity, it will stay 
there permanently. These special positions are called "libration points".

For the Earth-Moon mass ratio value of mu = 0.012155099, the x and y 
coordinates of these five equilibrium points are:

		x = +0.836892919    y = 0
		x = +1.155699520    y = 0
		x = -1.005064520    y = 0
		x = +0.487844901    y = +0.866025404
		x = +0.487844901    y = -0.866025404

Three of the libration points are on the x axis (collinear solutions) and the 
other two form equilateral triangles (equilateral solutions) with the 
positions of the Earth and Moon.

Program G3BODY will begin by displaying the following main menu: 

                Main Menu

	      < 1 > periodic orbit about L1

	      < 2 > periodic orbit about L2

	      < 3 > periodic orbit about L3

	      < 4 > user input of initial conditions

	      Selection (1, 2, 3 or 4)

The first three menu options will display typical periodic orbits about the 
respective libration point. The initial conditions for each of these orbits 
are as follows:

	(1) Periodic orbit about the L1 libration point

	X0 = 0.300000161        YDOT0 = -2.536145497
	MU = 0.012155092        TF = 5.349501906

	(2) Periodic orbit about the L2 libration point

	X0 = 2.840829343        YDOT0 = -2.747640074
	MU = 0.012155085        TF = 11.933318588

	(3) Periodic orbit about the L3 libration point

	X0 = -1.600000312       YDOT0 = 2.066174572
	MU = 0.012155092        TF = 6.303856312

Notice that each orbit is defined by an initial x component of position, X0, 
an initial y component of velocity, XDOT0, a value of Earth-Moon mass ratio 
MU, and an orbital period TF. The initial y component of position and x 
component of velocity for each these orbits is equal to zero. The algorithm of 
G3BODY is also valid for other three-body combinations.

If you would like to experiment with your own initial conditions, select 
option 4 from the main menu. The program will then prompt you for the initial 
x and y components of the position and velocity vector, and the value of the 
gravitational constant mu to use in the simulation. These requests are:
	
    Please input the x component of the radius vector
	
    Please input the y component of the radius vector
	
    Please input the x component of the velocity vector
	
    Please input the y component of the velocity vector
	   
    Please input the value for the Earth-Moon mass ratio

Here is a short list of initial conditions for several other direct and 
retrograde periodic orbits you may want to display with G3BODY.

	(1) Retrograde periodic orbit about m1 

		X0 = -.499999883        YDOT0 = 2.100046263
		MU = 0.012155092        TF = 11.99941766

	(2) Direct periodic orbit about m1 

		X0 = 0.952281734        YDOT0 = -.957747254
		MU = 0.012155092        TF = 6.450768946

	(3) Direct periodic orbit about m1 and m2 

		X0 = 3.147603117        YDOT0 = -.07676285
		MU = 0.012155092        TF = 12.567475674

	(4) Direct periodic orbit about m2 

		X0 = 1.399999991        YDOT0 = -9298385561
		MU = 0.012155092        TF = 13.775148738

The software will also ask for an initial and final time to run the 
simulation, and an integration step size. You might try 0.01 for the step 
size. The program will automatically adjust this step size to maintain the 
user-defined integration tolerance described below. These three prompts are:

    Please input the initial time
	 
    Please input the final time
	 
    Please input the integration step size
    (a value of .01 is recommended)

Please recall that these are "canonical" time units.

Main menu option 4 will also ask you to define the window for drawing the 
graphics. This window is specified by x and y coordinate pairs for the lower 
left hand corner and the upper right hand corner of the screen. These two 
prompts appear as follows:

    Please input the X and Y coordinates of the lower
    left hand corner of the graphics screen
	
    Please input the X and Y coordinates of the upper
    right hand corner of the graphics screen

These values are best determined by trial and error, and will depend on the 
actual orbit you are trying to plot. They are typically in the neighborhood of
+/- 1 to +/- 4. By changing the dimensions of the window, you can zoom in or 
out on the display screen. This allows you to see a better view of the motion 
as the spacecraft approaches the Earth or Moon, for example.

The program will also prompt you for a tolerance to use during the numerical 
integration. This number determines how accurately the trajectory is computed.
The request appears as follows:

    Please input the convergence tolerance
    (a value of 1D-8 is recommended)

The software begins by drawing the relative locations of the Earth, Moon and 
selected libration point on the graphics screen. Please note that the sizes of 
the Earth and Moon are not drawn to scale, and the actual physical dimension 
in the x and y directions will probably be distorted due to the aspect ratio 
of your monitor screen. 

During the graphics display, the distance between successive trajectory points 
and the program execution speed are indications of the integration step size. 
The step size should decrease as the satellite approaches the Earth or Moon. 
Can you explain why this happens?

