3D graphics in Resolver One using OpenGL and Tao, part II: an orrery

Posted on 17 September 2009 in 3D, Programming, Resolver One

In my last post about animated 3D graphics in Resolver One (the souped-up spreadsheet the company I work for makes), I showed a bouncing, spinning cube controlled by the numbers in a worksheet. Here's something more sophisticated: a 3D model of the planets in our solar system, also know as an orrery (click the image for video):

3D stuff:

To try this one out yourself:

A window showing the 3D animation of the solar system will appear when the file loads. To manipulate it, drag the mouse left/right to spin the solar system around the Y axis, or up/down to spin around the X axis (left to right in the initial state). Zoom in and out using the mouse wheel or +/-.

Here's how it works in more detail; I'll start with the spreadsheet and then go on to the OpenGL/Tao/IronPython support code that lives in OpenGLPlanets.py later.

The Planets worksheet

The Planets sheet of the spreadsheet is mostly static data, which (as I said in the video) is largely copied from Wikipedia. I shan't waste time here explaining exactly what an orbit's semimajor axis, inclination to the Sun's equator, longitude of the ascending node or argument of perihelion are because you can easily look them up -- or, even better, you can read this excellent tutorial on How to Calculate the Positions of the Planets, which I found halfway through building the spreadsheet and wish I'd found before I started.

The "Rings function" column in the Planets sheet uses the fact that Resolver One can store any kind of data in a cell to store one function per planet. Not the results of the function, but a pointer to the function itself, so that it can be passed around and called elsewhere by code that needs it. Just to make things a little bit more difficult for those who aren't comfortable with functional programming, the way I've implemented this is by writing a function that takes the parameters of the rings (internal and external radii and angle to the planet's orbit) and then returns a function that will draw rings with those parameters. So we have a cell that calls a function that returns a function that is then stored in that cell. All clear? Good.

(The only planet that has a rings function is, of course, Saturn. You might want to experiment with adding rings to other planets -- I think most of the gas giants have rings of some sort. In fact, perhaps a better name for this would be "Decorator function", because the function that is specified here can draw anything at all -- an interesting variation might be to modify it to draw the Earth's Moon. I'll leave that as an exercise for the reader...)

Moving briskly on to the other columns in the "Planets" sheet, we have the number of days in each "tick". When plotting an orbit, we divide it into a given number of points (700 is the default); a "tick" is the amount of time it takes the planet to get from one tick to the next. (The first version of this orrery plotted a point for every (Earth) day, so there were 88 ticks in Mercury's orbit, 224 or so in Venus', 365-odd for Earth, and... well, things got a little hairy when it came to Pluto's 248-year orbit.) Anyhow, this number is used when calculating the orbits. An interesting side-point -- instead of being calculated by formulae, it's calculated by a loop in the pre-formulae user code. This is to handle a rather complex dependency issue that I shan't explain here, though if anyone's interested I can give details in the comments.

After the "days/tick" column, we just have a few columns whose sole purpose is to convert the user-friendly degrees in which we have entered some of the orbital parameters into math-function-friendly radians. These are done in the same pre-formulae code as the "days/tick" calculations, for the same dependency-related reason.

The other worksheets

Now let's move on to the worksheets for the planets. They're all basically the same, so let's look at the one for Mercury as an example. Working through the columns one by one:

And that's it for the worksheets. Now, on to the user code.

The pre-constants user code

Look at the blue section in the code editor that is below the spreadsheet grid for this code. It starts off with a definition of the number of points we want to calculate per orbit. The default of 700 means that it takes 14 seconds to calculate everything for every planet on my (not very speedy) machine; you can increase this to get smoother animations with slower calculation times, and decrease it if you don't mind the planets leaping their way around their orbits.

The code then goes on to define the functions that do the heavy lifting for the column-level formulae that we looked at earlier in the per-planet worksheets. calculateMeanAnomaly, calculateTrueAnomaly, calculateOrbitalDistance, calculateX, calculateY, and calculateZ are all simple enough (and did I mention that there's a great web page describing how they work?), but calculateEccentricAnomaly is a little complicated. The problem is that we have a value called the mean anomaly (M), another called the eccentricity (ε), and we want to work out the eccentric anomaly (E), but the equation that relates them is:

M = E - εsin(E)

If you can't see how to rearrange that so that it tells you how to get E, given a value for M, you're not alone. The equation is transcendental, and there is no simple solution. Instead, we have to use a numerical method to work out a best-guess value for E. calculateEccentricAnomaly uses the Newton–Raphson method to do this, generating a result of E that is correct to within 5 decimal places.

The pre-constants code finishes off with the function used for rings; it takes the angle to the planet's orbit, and inner and outer radii for the rings, and returns a function that draws the rings. This function in turn takes a parameter that tells it how much to scale the radii -- this value is calculated by the main OpenGL code in OpenGLPlanets.py, which we'll come to later.

The pre-formulae user code

Further down the code editor window you will find this user code section, with a green background. It calculates various values that the formulae in the spreadsheet will use.

It consists of a loop through the worksheets defining individual planets; for each, it:

Once this is all done, there is enough information for the column-level formulae in the planets' worksheets to calculate the detailed orbits.

The post-formulae user code

This is very similar in structure to the code in my last post; we import the OpenGL code used to display the planets, check if we already have an open graphics window, and create one if we don't.

Next, we clear the window's list of planets to draw, and then iterate over the worksheets that contain planet data. For each one we create a new Planet object populated with the calculated positions for the orbit and with the other details of the planet -- radius, colour, ring function, and so on -- and add that to the window's list of planets. Once that's done, the window can display the orbits.

That's everything in the spreadsheet. Now, let's look at the supporting OpenGL code in OpenGLPlanets.py.

The IronPython support code

I won't go into this is huge amounts of detail; there are many good OpenGL tutorials and references out there, and duplicating stuff that's already been done much better would be a waste! In addition, much of what's in there is the same as the code in my last post, and was explained there. So, a quick run through the new bits in the code:

And that's it! Hopefully that was all clear, but if anything needs more explanation, please do leave a comment and I'll clarify things as much as I can.