Orbital Mechanics Introduction

• Download demo project - 12.7 MB

Introduction

This primary project is called OrbitMech. Two other Visual Studio projects are described in detail below. OrbitMech is a Windows Forms project written in C# using Windows Forms and OpenGL. It allows you to change orbital parameters to visualize the orbit of a small body around a much larger body, such as the orbit of a satellite around the Earth. The parameters you can modify are described below and are:

1. Semi-major Axis
2. Eccentricity
3. Inclination
4. Longitude of the Ascending Node (Omega)

Background

The subject of orbital mechanics is very complex, but a simple subset known as the "two body problem," in which an object is in orbit around an object that is significantly more massive involves fairly simple mathematics. A satellite orbiting the earth, a spacecraft orbiting the moon, or a planet orbiting the sun are examples. That is, we assume that the mass of the smaller object is negligible compared to the mass of the large object. The First Law of Planetary Motion, discovered by Kepler, states that the shape of an orbit is an ellipse with the massive body at one focus of the ellipse. Newton later developed the mathematics of elliptical orbits. The more massive body will be referred to as the "planet."

An ellipse, as you will recall from your geometry class, is the set of planar points such that the sum of the distances (shown as dashed lines) from two fixed points, known as foci (singular: focus) is constant. An ellipse is defined by the Semi-major Axis and the Eccentricity. The major axis of an ellipse is its longest diameter and is a line segment that runs through the center and both foci; the Semi-major Axis is half the major axis. The eccentricity has a value from 0 to 1 and is a measure of how "non-circular" the ellipse is; a circle is an ellipse with an eccentricity of zero. In the image to the left, the ellipse has an eccentricity of about 0.6. The shorter axis, the minor Axis, is also shown.

There are six parameters that completely define an orbit, four of which will be addressed in this project. The first two orbital parameters define the shape of an ellipse, and therefore the shape of an orbit. They are the Semi-major Axis and the Eccentricity of the ellipse.

In the image above, the planet (the massive object, Earth in this example) is shown along with the orbit (in blue) in a Windows Form. The XYZ axes are aligned with zero degrees longitude, ninety degrees longitude, and the North Pole, respectively. At the bottom of the Form are four track bars. The values of the Semi-major Axis and the Eccentricity can be changed by using the corresponding Track Bars. Note that as you change the eccentricity, the semi-minor axis of the ellipse is displayed in parentheses next to the value of the semi-major axis. The plane that contains the orbit, the Orbital plane, is shown in gray. The Major Axis is shown in cyan and the Minor Axis in magenta.

The next two orbital parameters define the orientation of the orbital plane. They are the Inclination, which is the vertical tilt of the ellipse with respect to the orbital plane measured at the ascending node (where the orbit passes upward through the orbital plane,) and the Longitude of the Ascending Node, which horizontally orients the ascending node of the ellipse. The Longitude of the Ascending node is commonly referred to by the Greek letter omega, Ω. These two parameters can be changed by using the corresponding Track Bars. (For simplicity, I've chosen to not use the final two orbital parameters: the Argument of periapsis and Mean anomaly.) The orbital period (in minutes,) apogee (highest point in the orbit) and altitude at apogee are displayed in the lower left side of the form.

Using "Display" on the menu bar, you can toggle the drawing of the axes, planet, and orbital plane. The "camera" can be moved by using the left, right, up, down, plus and minus keys:

• Up/Down Arrow: tilt camera by phi (F) degrees
• Right/Left Arrow: rotate camera by theta (?) degrees
• Plus/Minus: zoom camera In/Out by R units

The three Camera settings are displayed in the lower right side of the form. The calculation of phi, theta and R is contained in the Camera class which uses the OpenGL utility method gluLookAt().

OrbitMech Configuration File

The 4 orbital parameters (Eccentricity, Inclination, Semi-major Axis, and Longitude of the Ascending Node or Omega,) as well as the camera settings can be specified in the App.config file. For example, here is a snippet of an App.config to specify an Eccentricity of 0.05, an Inclination of 50 degrees, a Semi-major Axis of 6556 kilometers, and an Omega of 30 degrees:

<setting name="Eccentricity" serializeAs="String">
	<value>0.05</value>
  </setting>
<setting name="Inclination" serializeAs="String">
	<value>50</value>
</setting>
<setting name="SemiMajorAxis" serializeAs="String">
	<value>6556</value>
 </setting>
<setting name="Omega" serializeAs="String">
	<value>30</value>
   </setting>

Here is a snippet of App.config to specify a camera tilt (phi, F) of 20 degrees, a rotation (theta, ?) of 60, and a zoom of 27:

<setting name="CameraPhi" serializeAs="String">
	<value>20</value>
</setting>
<setting name="CameraTheta" serializeAs="String">
	<value>60</value>
</setting>
<setting name="CameraR" serializeAs="String">
	<value>27</value>
</setting>

The background color and color of the orbital plane can also be specified in the App.config file, as the ARGB values or color name, for example:

<setting name="OrbitPlaneColor" serializeAs="String">
	<value>26, 224, 224, 224</value>
</setting>
<setting name="BackgroundColor" serializeAs="String">
	<value>Gray</value>
</setting>

When you Exit OrbitMech, you will be asked if you want to save your orbit parameters and camera settings in the Applications Settings (in AppData\Local on Windows 7,) and on starting OrbitMech, you will be asked if you want to restore those settings. See SaveValues() and RetrieveValues() in OrbitMechMainForm.cs for details. Alternatively, the orbital parameters and camera settings can be saved in a text file using "File->Save As..." from the menu bar, and restored by using "File->Open...".

Setting the Planet

A key constant in computing the position and speed is the product of the Gravitational constant and the mass of the planet, the Standard Gravitational Parameter, commonly referred to by the Greek letter Mu, μ. You can change this value using the configuration item Mu in the App.config file. Mars, for example, has μ = 42828. You can also change the texture that is applied to the sphere using a bitmap stored in the Textures directory. The Textures directory and contents must be copied to the bin\Debug and bin\Release directories. The planet texture is set with the configuration item PlanetTexture in the App.config file. For example, setting it to "mars.jpg" will paint the sphere with a bitmap of Mars. (Thanks to http://vasilydev.blogspot.com for the code that paints the bitmap on the sphere.) You can set the radius (in kilometers) of the planet about which the body is orbiting with the configuration item PlanetRadiusKM in the App.config file.

Here is a snippet from App.config showing how to set Earth as the planet:

<setting name="PlanetTexture" serializeAs="String">
	<value>earth.jpg</value>
</setting>
<setting name="PlanetRadiusKM" serializeAs="String">
	<value>6371.1</value>
</setting>
<setting name="Mu" serializeAs="String">
	<value>398600.0</value>
</setting>

Here is a snippet from App.config showing how to set Mars as the planet:

<setting name="PlanetTexture" serializeAs="String">
	<value>mars.jpg</value>
</setting>
<setting name="PlanetRadiusKM" serializeAs="String">
	<value>3389.1</value>
</setting>
<setting name="Mu" serializeAs="String">
	<value>42828.0</value>
</setting>

Here is a snippet from App.config showing how to set the Moon as the "planet":

<setting name="PlanetTexture" serializeAs="String">
	<value>moon.jpg</value>
</setting>
<setting name="PlanetRadiusKM" serializeAs="String">
	<value>1737</value>
</setting>
<setting name="Mu" serializeAs="String">
	<value>4902.8</value>
</setting>

Setting Altitudes

If you would like to enter the altitudes of the orbit and have OrbitMech compute the Eccentricity and Semi-Major Axis, set both the Eccentricity and Semi-Major Axis to -1 (minus one) and specify the altitudes in kilometers as in this snippet from App.config:

<setting name="Eccentricity" serializeAs="String">
	<value>-1</value>
</setting>
<setting name="SemiMajorAxis" serializeAs="String">
	<value>-1</value>
</setting>
<setting name="MaxAltitude" serializeAs="String">
	<value>297</value>
</setting>
<setting name="MinAltitude" serializeAs="String">
	<value>111</value>
</setting>

Apollo 8, the first manned lunar mission, had an initial orbit with an inclination of 32 degrees and the minimum and maximum altitudes were 297 and 111 km. Using those values as shown in the above App.config snippet, OrbitMech computes a semi-major axis of 1941 km and an eccentricity of 0.048 as shown in the image below. (Note how I have turned off the Axes using "Display->Axes Off", and set the Background and Orbit Plane colors in App.config.)

Using the Code

Download and extract the project, and open OrbitMech.sln with Visual Studio. It consists of three projects:

1. OrbitMech - A Windows Forms project to display an orbit in 3-D using OpenGL
2. KeplerForm - A Windows Forms project to animate an object in orbit around a planet in 3-D using OpenGL. It uses KeplerLib to perform the calculations.
3. KeplerLib - A library to numerically solve Kepler's equation (For an introduction to numerical methods, see My numerical methods project.)

Build the solution using "Build"->Rebuild Solution" from the Visual Studio Main Menu. Use Ctrl-F5 to Start without Debugging and display a Windows Form as shown in the image above.

OpenGL

Although this article is not intended as an OpenGL Tutorial, some explanation is provided here for those who are not familiar with it. The OpenGL Programming Guide, or "Red Book" as it is known, is a great resource for learning OpenGL.

As an example, here is the code from DrawPlane() in OrbitalBody.cs to draw the orbital plane. With OpenGL, all geometric objects are ultimately described as an ordered set of vertices bracketed by glBegin() and glEnd() (note how I've used C#'s curly braces so you can match glBegin() and glEnd() calls easily with Visual Studio's CTRL-] - this is simply a convenience and is not required by OpenGL). In this example, the orbital plane is a rectangle with four vertices. The dimensions of the rectangle are defined by the semi-major and semi-minor axes of the ellipse. The rectangle's corners are defined by the four glVertex3f() calls and the rectangle lies in the X-Z plane since the Y coordinate (the second argument) in the call to glVertex3f() is zero.

Gl.glBegin(Gl.GL_QUADS);  // Each set of 4 vertices form a quad
{
...
Gl.glVertex3f(-semiMajorAxis + focus, 0.0f, -semiMinorAxis);
Gl.glVertex3f(+semiMajorAxis + focus, 0.0f, -semiMinorAxis);
Gl.glVertex3f(+semiMajorAxis + focus, 0.0f, +semiMinorAxis);
Gl.glVertex3f(-semiMajorAxis + focus, 0.0f, +semiMinorAxis);
}
Gl.glEnd();

A key concept of OpenGL is that of the modelview matrix to apply a Modeling Transformation, for example rotating or translating part of a scene. As an example, here is the code from DrawEllipse() in OrbitalBody.cs to draw the actual orbit. When the inclination of the orbital plane is changed by moving the Inclination Slider, the value (in degrees) of the inclination obtained from the slider control is used as the first argument to glRotatef(), a Modeling Transformation method. The next three arguments specify the x, y, and z components of a vector about which the rotation occurs, in this case the X-axis. The code to rotate by omega is a similar call to glRotatef().

Gl.glPushMatrix();
{
Gl.glRotatef(orbitalParameters.Inclination, 1.0f, 0.0f, 0.0f);
		...
DrawPlane(semiMajorAxis, semiMinorAxis, (float)focus);
	...
}
Gl.glPopMatrix();

As the Red Book says:

In effect, glPushMatrix() means "remember where you are" and glPopMatrix() means "go back to where you were."

So the above code has the effect of rotating the orbital plane - that is, changing its inclination - without rotating anything else in the scene. Again, I've used C#'s curly braces so you can match glPushMatrix() and glPopMatrix() calls with CTRL-] and this is simply a convenience and is not required by OpenGL.

Tao Framework

OrbitMech uses the Tao Framework OpenGL implementation. These files are included in the OrbitMech project and must be copied to the bin\Debug and bin\Release directories:

1. Tao.OpenGL.dll
2. ShadowEngine.dll
3. glut32.dll

Numerical Solution to Kepler's Equation

In order to compute the position and speed of the body as it orbits, it is necessary to numerically solve Kepler's equation.
I am grateful to Richard Gonsalves for his implementation of a Runge-Kutta numerical integration of Kepler's equation which I've converted from C++ to C# and placed in a library called KeplerLib. I modified it to iterate for one orbit, and to display distances in km and velocities in km/sec. Right-Click on "Kepler" in the Visual Studio Solution Explorer, then click "Set as Startup Project." Use Ctrl-F5 to Start without Debugging and display a Windows Form as shown in the image above. The satellite's orbit around the Earth has the following orbital parameters:

• Eccentricity = 0.6
• Semi-major axis = 18125 km
• Inclination = 0
• Omega = 0

(In this image, I've used the down arrow key to set the Camera's Phi (F) value to 65 degrees. The Camera's Phi, Theta, and R values are displayed in the console window.)

Clicking the "Start" button will start an animation illustrating a satellite orbiting a planet. (After pressing "Start", the button label will change to "Pause". I pressed "Pause" in order to get the screen shot. Pressing "Continue" will continue the animation.) The constructor KeplerForm calls the library method KeplerCalc.Integrate() which in turn solves the equations of motion and creates a list of state vectors, List<StateVector> svList. Each state vector represents the position and velocity of the orbiting body at an instant of time. This list of state vectors with their position and velocity of the satellite is used by the KeplerForm's DrawOrbit() method. As it iterates through each StateVector in svList, DrawOrbit() draws the satellite (a small red square) at the StateVector's position, a PointF struct. A timer called tmrPaint is used to drive the animation. The timer's interval is computed based on the speed of the satellite, so in the animation the satellite will speed up as it reaches the closest point in its orbit to Earth (perigee), and slows as it reaches the furthest point in its orbit (apogee,) a phenomenon observed by Kepler himself.

Console Window

As with my NXT Bluetooth Monitor project, I've attached a console to the Windows Forms projects for debugging purposes using AllocConsole and .NET Interop. It is enabled by default; to disable it, simply comment out the first 7 lines of Main() in OrbitMechMain.cs or KeplerFormMain.cs:

IntPtr ptr = GetForegroundWindow();
int pid;
GetWindowThreadProcessId(ptr, out pid);
Process process = Process.GetProcessById(pid);
AllocConsole();
AttachConsole(process.Id);
Console.WriteLine("{0} OrbitMech started", DateTime.Now);

History

• Version 2.0.0.0 Omega is now shown as East or West longitude (Click Help to see the version)