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