Last edited May 4, 2008 by Will Davis
Find this document at http://cheer.math.edu/reulab/web/wldavis
illiSol: A Modern Orrery - 2008
Abstract
     The goal of this project is to create a modern "orrery"
- that is, a model of our solar system - for immersive virtual environments
(specifically, the CUBE and the CAVE here at the University of Illinois), using
OpenGL and the virtual environment development tool Aszgard. As of April 3,
2008, the project is expected to have 3 base modes, where a base mode
designates how the planetary system is viewed.
         Base Modes:
         Base Mode 1: This
viewing mode will be designated as the initial mode. It should feature a
moderately accurate and scaled view of the solar system. After primary
functionality is available, further secondary functionality may be added for the
wand (virtual environment control device) to adjust the scale of the planets to
be more visible. This viewing mode should also feature the option to have the
orbits of each planet toggled.
         Base Mode 2: This
viewing mode is, in reality, more of a "flying mode". In this mode, the user is
given full control, via wand and headtracking, of their direction. The
specifics of this mode are still being determined, although it is likely to
posess the interesting feature of variable time - that is, time will most
likely vary according to proximity or some power of the proximity to objects of
interest (which will be predetermined). This is due to the impractical time it
would take to actually travel through our solar system in linear time. By final
presentation, the user should be able to highlight a destination (to aid them
in orienting themselves toward it), and be able to see it, even when it is very
far away.
         Base Mode 3:
(Developmental) This viewing mode may not come to fruition, as it is a nice
touch but may not be feasable in the time constraints. Because Jupiter and its
moons are often of interest to astronomers, this viewing mode would present an
explorable region bounded around Jupiter, but enough to fly from moon to moon.
     This project is being developed in concert with Mike
Hutches' project, a program to simulate the Voyager Mission through the Sol
System. This addition may add the opportunity for an Auxiliary Mode, which is
similar to Base Mode 2 but follows a predetermined path - namely, the path of
the Voyager - with the same feature of variable time present in the system.
Also, the data collected under Mike's project will be essential in initializing
the solar system accurately to a given time, which, as a further option, may
be provided by the user (although it will most likely be restricted to the
past). Reasonably, it should follow that the planets in each viewing mode will
move in their correct orbits according to time.
Here are some examples of pictures that will be mapped to their planetary
spheres:
       Earth
       Jupiter
       Saturn
Update 4-28-08: It turns out that for our image mapping to work onto a
gluSphere() (in its simplest form), we require .raw image files, not .png or
.jpg (although this is possible with extra libraries). Thankfully, 1
UIC has a tutorial/homework assignment for just such a thing with all the
necessary .raw files, and a description of how to use them.
Update 5-4-08: It turns out this method of texture mapping is very
inefficient, as it creates a set of mipmaps for every texture for every step the
texture is drawn, which results in quite a slowdown of the program. This is
currently the greatest barrier in advancing illiSol to the "next level." If any
more work is to be done, it should be to integrate the texture mapping with
the Syzygy environment so that it can be done *far* more efficiently, and pave
the way for more improvements on the features and structure of the program.
Some Mathematics Involved
Orbital Elements Introduction
     Fortunately for us, a precise guide to calculating
planetary positions is available online, courtesy of Jet Propulsion Lab at
CalTech2. The document provided has a step-by-step guide to
calculating the necessary orbital elements, which, for the most part, are
referenced similarly in the Wikipedia article linked to above. Following that,
it is possible to calculate each planet's position a couple different useful
coordinate systems. The coordinate system of interest to us, of course, is the
heliocentric system, a three-dimensional Euclidean system with a central body
(namely Sol) at the origin. So what exactly is it that we are doing?
     It is fairly commonly known that planets generally follow
a slightly elliptical orbit around the central body. An exception is Mercury,
which has a much higher orbital eccentricity, and the comets in our system,
which have a very high eccentricity. However, the plane that these bodies orbit
in are not the same for each orbiting body - which is why we need a third
coordinate to describe its "inclined" position. However, dealing with
non-uniform orbits as such introduces a certain amount of complexity. The way
we can deal with this complexity is solving it in a piecewise fashion.
     By now we have determined that the plane of orbit is
inclined from a reference plane that the central body lies in. For
simplicity's sake, we will think of this reference plane as the xy plane
in the xyz coordinate system. The plane of orbit then, must contain both
the central body and the orbiting body. I will refer to this as the orbital
plane, while the angle separating the orbital plane and the reference
plane is commonly known as the angle of inclination.
     Because each orbit is unique and elliptical, we want to
orient it in a useful way. To do this, we need to know a little bit about the
characteristics of elliptical orbits (this next section assumes some prior
knowledge of ellipses).
Facts:
     1) It has turned out that the central body will
always lie on one of the foci of the ellipse.
     2) Given this, the points of greatest and least
distance from the central body to the path of the orbit will always
lie on the same line segment. The point of greatest distance, when considering
heliocentric orbits, is known as the aphelion while the point of least
distance in the same context is known as the perihelion.
     We can see, just from visual observation of our solar
system, that not all orbits are oriented in the same way. Often, we define a
sort of longitude in the reference plane, and we need a zero longitude. In our
solar system, and hence in this project, the zero longitude in the reference
plane is given by the point of vernal equinox, the point in the Earth's orbit
when the equatorial plane of the Earth is in orbital plane and the Sun is
shining directly at equator. For satellites of orbiting bodies, the reference
direction is, for general purposes, parallel to the Vernal equinox vector.
Interestingly, JPL has given us what they call the
longitude of perihelion, which, I believe, given the Wikipedia diagram,
is the angle that would be obtained if you begin at the reference direction
vector and travel counter-clockwise along the reference plane until the
the ascending node is reached, at which point the angle then continues up
the inclined orbital plane, still counter-clockwise, to the perihelion vector.
JPL also provides the longitude of ascending node, which spans from the
reference direction vector to the ascending node, in a counter-clockwise
fashion. This way, the argument of perihelion can be obtained by
subtracting the longitude of ascending node from the the longitude
of perihelion:
     The next two orbital elements used in position calculation
deal with a sort of approximation. These are the mean longitude and
the mean anomaly, which are calculated as if the orbit were actually
circular and the inclination angle were zero. First, the mean anomaly is the
analogue of the true anomaly under such ideal circumstances. It can be
calculated, with some error for certain ranges of dates, by subtracting the
longitude of perihelion from the mean longitude. The mean longitude, which has
been supplied in the JPL documents, is the longitude of the orbiting body in
the plane under ideal circumstances.
     To continue the calculation, we need to find the
eccentric anomaly, which requires an iterated approximation method to
compute, as shown in the JPL documents. The eccentric anomaly represents the
angle obtained in the following manner:
     1) Circumscribe the orbital plane with a circle in
the same plane.
     2) In the orbital plane, draw a line from the
planet's position down to the semi-major axis of the orbit so that it is
perpendicular to the semi-major axis.
     3) Extend that perpendicular line from the planet
away from the semi-major axis to the point where it intersects the auxiliary
circle that we have circumscribed the ellipse with.
     The resulting angle from the semi-major axis to the point
where the line intersects the auxiliary circle is the eccentricity angle E
.
Wolfram's
Mathworld provides a good visual aid and explanation for this, where point
P is the planet in question. The JPL documents approximate this using iterated
computational methods.
     Now that we have the eccentric anomaly, we can finally
start to compute our coordinates. First, we find the heliocentric coordinates
in the orbital plane, with the central body at the origin and the x'
axis aligned towards the perihelion. This can be accomplished using a little
trig, the eccentricity and semi-major axis of the orbit, and the eccentric
anomaly we just computed. Because these are simply planar coordinates, the
z' coordinate will be zero at all points.
     We can then translate this position the
reference plane, which the JPL documents refer to as the J2000 ecliptic plane.
But once again, this is fairly straightforward since we know all the relevant
angles describing our orbit. These equations will translate our planetary
position for us:
     See the JPL Documents (linked below) for more equations
and detail.
     And that's it! We now have compact way to calculate the
position of any of the planets in the Sol System, and even some of their
satellites, given only the date in a range from 3000 B.C. to 3000 A.D.!
Mode of Operation - OpenGL
     It just so happens that the way this application is
structured seems to mesh perfectly with how OpenGL works. To demonstrate this,
let's provide a typical example. Let's say we have a planet that is in orbit
around a central, stationary (for the purposes of this application) star. For
further interest, let's say that our planet has a moon, as in the crudely-drawn
diagram above. Let's examine how our system moves, starting with the smallest
element - the moon.
     1) The moon may move in two ways - it may rotate
around its own central axis and it is also in orbit around the planet.
     2) The planet also moves in two ways - it rotates
around its own central axis and it is also in its own orbit around the
stationary star.
     However, as one can see, because of 2), the moon is
also moving around the star, in much more complicated path than any
elliptical orbit. But no worries, we can still have satellites in this
application, and here's why.
     The way OpenGL works is through, essentially, matrix
multiplication. The OpenGL design takes each pixel, applies a set of matrices
to it, and immediately sends it out to the display (see Dr. George Francis'
paper on Real-time Interactive Computer Animations
here
for deeper details)3.
     Essentially then, we theoretically can have infinitely
many levels of orbit. In other words, we can have satellites that have
satellites that have satellites that have satellites... and so on, where the
only limit would be the amount of stack memory allocated to the OpenGL process.
This is because we apply each transformation matrix one at a time. For
instance, we apply a single transformation matrix for a planet to put it in its
orbit, and we then draw the planet. After that, we apply another transformation
matrix and then draw the satellite. In this case, the planet would only have
one transformation matrix applied to it, while the satellite would have both
matrices applied to it. This is possible because of the stack nature
of OpenGL, a Last In, First Out (LIFO) data structure.
     As seen in the Mathematics section, we first obtain the
heliocentric coordinates of each planet it its own orbital plane before
we get the heliocentric coordinates in the reference plane. So given our
Mode of Operation, is obtaining the heliocentric coordinates in the reference
plane even necessary? No, it's not! The design of OpenGL allows us to eliminate
a level of computation from our calculations while at the same time accounting
for each successive transformation. Indeed, this should be a very useful trick!
Progress Reports
Friday, April 18, 2008 - The initial programming of illiSol has begun, as of
Wednesday, April 16, 2008. I have created a header file for the program, KeplerianElements.h, that holds all the data we might want
to compute the positions of the 8 major planets (since Pluto is no longer part of the
club). The main development goal right now is to get the Orrery working in OpenGL - that
is, I want the system to "be there" before I try to do anything else. The current
illiSol.cpp file is mostly a skeleton, but the main mathematics for compuing
planetary positions are present.
Monday, April 28, 2008 - In the past week and a half, many changes have been made to the core structure of illiSol. I abandoned my endeavor to create
my own full-fledged OpenGL/Syzygy application from scratch, opting instead to
use the skeleton file included with the Aszgard development kit, which is basically a shell program with lots of useful OpenGL and Syzygy tidbits. As of today,
I now have a fully functional standalone copy of illiSol that works on Windows XP and Windows Vista, but no attempts have been made to compile it for *nix. Fully functional includes the following:
     - All of the major planets in their proper orbits (sorry Pluto, you were the weakest link).
     - A system that changes with time - the rate of change of
time is variable, but always increasing. The initial date for now is August 20, 1977, around the time of the Voyager 2 Launch.
     - All of the planets, and the Sun, now have surface images mapped to their respective spheres. Saturn also has a few rings to distinguish it.
     - I now also have the ability to toggle magnification on and off. The planets are magnified by 100x and the sun is magnified by 5x. The distances between the planets remains the same, regardless of the state of magnification. Any larger value of magnification for Sol would engulf Mercury.
     - The most recent touch was to add the date to the window text. The user can now see the Julian Day Number as well as the current date in the Gregorian Calendar, to get a sense of how quickly time is advancing.
     I have also decided to remove the option of drawing the orbits of the planets into the Solar System. This becomes fairly memory intensive and and at least doubles the computational load on the machine. Since the application itself is already fairly resource intensive and not running as quickly as I would like, I have decided to focus my final efforts on optimizing the current computations and graphics rather than adding to the problem, which will most likely be magnified in the CUBE and CAVE anyhow. Also, because of time constraints, it will not be possible to complete viewing mode 3 (reference above).
     Here are some screenshots of the progress thus far this was taken running on a quad-core processor with 4GB of RAM, running Windows Vista. You can see the frame rate as around 20, whereas, for comparison, on a normal dual-core processor with 2GB of RAM running Windows XP, the frame rate ranges from 7 to 12 fps.
illiSol: Earth
illiSol: Magnification Off
illiSol: Magnification On
Wednesday, April 30, 2008 - There is now a Windows executable (with data
files) and source code available. NOTE: To run the executable you will need
the required .dll files, which are not included in the package. A standalone
with the required .dll files may be made available in the future. Regardless,
if you have Aszgard, you should have the files you need, and there are
instructions for compiling included in the source code zip.
To download, right click on the link and choose "Save Link As...":
Windows Executable and Data Files
Source Code, Compiling Instructions, and Data Files
Friday, May 9, 2008 - New updated executables and source code are
available! This version features improved texture mapping efficiency and
correctly mapped buttons for the virtual environment controllers.
To download, right click on the link and choose "Save Link As...":
Windows Executable, Data Files,
and .dll files - illiSol can now be run on any Windows machine without
Aszgard.
Source Code, Compiling Instructions, and Data Files
Reference Materials
1"Texture Mapping... A Quick Overview." Andy's shell program.
http://www.evl.uiuc.edu/julian/cs488/2005-11-10/index.html
2Standish, E.M. "Keplerian Elements for Approximate Positions of
the Major Planets." Jet Propulsion Laboratory, California Institute of
Technology.
http://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf
3Francis, George K., Ph.D. "Metarealistic Renderings of Real-time
Interactive Computer Animations." University of Illinois, Nov. 7, 2001.
http://www.math.uiuc.edu/~gfrancis/public/mrtica13may.pdf