Object Oriented Astronomy–Part 1

Sunday January 26 2014 - , , , , , , , , ,


One of the things that bugs me about astrometry and astrophysics libraries is that, more often than not, they are direct ports from some original Fortran 77 version. The ported code tends to be very ‘procedural’, which forces the user of the library into a certain one-dimensional way of thinking. You may think there’s nothing wrong with that, and you may very well be right, but I think we can do better. I think its about time that someone stepped back from the spaghetti code and did a proper object-oriented analysis and brought astronomy programming out of the 1970s. This is a series of blog posts I’ll be publishing to cover my attempt at doing so.

Where Do We Start?

In this first instalment, I’ll cover some background information and try to explain why I’m doing this. I wanted to start out with something useful, like calculating the position of a planet for a given date and time. I want to work test-first (using MSpec), try to be as SOLID as I know how, and produce clean code.

I’ve examined and been inspired by several libraries, such as

  •  NOVAS (the US Naval Observatory Vector Astrometry Software, originally written in Fortran and ported to C and Python)
  • The ASCOM rendition of the above, which is a Visual Basic (.net) port
  • AA+, a C++ class framework for computational astronomy by P. J. Naughter which in turn is based on algorithms presented in the book Astronomical Algorithms by Jean Meeus.
  • AAPlus, a C# adaptation of AA+ by ‘ezz’.

Of the above, AAPlus has probably been the most instrumental in prodding me into action. The code is at least in C# (my target language) and has a unit test suite that compares its results against AA+, which all bodes well. However, it provides little advantage over the original AA+ other than simply using different language syntax, it would almost have been easier just to write a NativeMethods class full of DllImport declarations.

AAPlus uses something called VSOP87 to compute orbital positions, as does AA+. NOVAS appears to use an interpolation technique based on JPL planetary ephemeris data files. While it makes no pretense at being object oriented, never progressing beyond C, it is particularly well documented. ASCOM relies on NOVAS-C and also has something called the ‘Kepler Orbit Engine’ which can compute positions of minor planets and comets as well as planets. The idea of an ‘engine’ for computing orbital positions is something that resonates with me and I’ll probably carry this through into my final solution.

One very useful resource I discovered is a web site that will actually generate VSOP87 code for a given planet. This is helpful for two reasons.

  1. It is a useful reference implementation for validating results
  2. It provides an excellent example of what I think is wrong with many of the existing libraries.

What’s Wrong With Existing Implementations?

I don’t really see it as a matter of right and wrong. All of the solutions hitherto mentioned serve a well defined purpose and work. In that sense, there is nothing ‘wrong’ with them other than I don’t like the way things are structured. Perhaps that is more to do with me than these implementations, but indulge me for a while I arbitrarily and unfairly pick on the code generator web page as an example of what I don’t like. I used it to generate VB.net code as that is the closest to my target of C# (conversion to C# is trivial using the Telerik Code Converter).

Before going further, I need to explain a bit about VSOP87 and how it works. It is essentially an analytical technique for computing orbital positions. I will not go too deep into the theory, It is complicated to understand but fortunately complicated maths often boils down to a relatively simple algorithm. If you are interested in the mathematical background, try Numerical expressions for precession formulae and mean elements for the Moon and the planets and for implementation details try http://www.astrosurf.com/jephem/astro/ephemeris/et300VSOP_en.htm and for a general overview see http://en.wikipedia.org/wiki/VSOP_(planets).

Cutting to the chase, computing an orbital position with VSOP87 involves calculating the value of a formula of this ilk:

X, Y, Z are (in this version) the Cartesian coordinates to be determined; Ai, Bi, Ci are essentially constants read from records in a data file, and T is in thousands of Julian days from JD2000, T = (jd - JD2000) / 365250. The exact formula depends on the results required. VSOP87 has 6 different variants, which generate:

  • VSOP87 - Heliocentric ecliptic orbital elements for the equinox J2000.0; the 6 orbital elements, ideal to get an idea how the orbits are changing over time
  • VSOP87A - Heliocentric ecliptic rectangular coordinates for the equinox J2000.0; the most useful when converting to geocentric positions and later plot the position on a star chart
  • VSOP87B - Heliocentric ecliptic spherical coordinates for the equinox J2000.0
  • VSOP87C - Heliocentric ecliptic rectangular coordinates for the equinox of the day; the most useful when converting to geocentric positions and later compute e.g. rise/set/culmination times, or the altitude and azimuth relative to your local horizon
  • VSOP87D - Heliocentric ecliptic spherical coordinates for the equinox of the day
  • VSOP87E - Barycentric ecliptic rectangular coordinates for the equinox J2000.0, relative to the barycentre of the solar system.

I’ve chosen to use VSOP87B for this exercise, which generates spherical coordinates, i.e. latitude L and longitude B (in radians) and radius R (distance in Astronomical Units, AU). Note the two sigmas in the formula. There are 5 main terms each of which is itself a series of n terms which are affected by the value of T from the ‘parent’ term. All of that has to be performed three times, once for each component of the final coordinate. This ends up chomping through a LOT of data. For Earth, the latitude  L(α=0) term has 622 rows of data; L(α=1) has about 377 rows, L(α=2)  has 143 rows, L(α=3) has 22 rows, L(α=4) has 11 rows and L(α=5) has 4 rows, a total of some 1,179 rows of data that needs to be read, computed and summed. That is repeated (with different data) for longitude (B) and Radius (R), so upwards of 3,000 rows of data need to be processed. The code generator encodes that data like this:

   1: Public Function Earth_L0(ByVal t As Double) As Double ' 623 terms of order 0
   3:     Dim L0 As Double
   4:     L0 += 1.75347045673
   5:     L0 += 0.03341656453 * Math.Cos(4.66925680415 + 6283.0758499914 * t)
   6:     L0 += 0.00034894275 * Math.Cos(4.62610242189 + 12566.1516999828 * t)
   7:     L0 += 0.00003417572 * Math.Cos(2.82886579754 + 3.523118349 * t)
   8:     L0 += 0.00003497056 * Math.Cos(2.74411783405 + 5753.3848848968 * t)
   9:     L0 += 0.00003135899 * Math.Cos(3.62767041756 + 77713.7714681205 * t)
  10:     ' Continues up to 622 rows (elided to avoid terminal boredom)

Yuck! First, it’s clear that what this generator is really doing is generating _data_, not code. As the web site states:

In practice, writing the full-precision, hard-coded VSOP87 source code manually would be too prohibitive a task for a single human to attempt, due to the multiple thousands of mathematical terms, and very easily prone to typographical errors

Very true, but the original, unadulterated data files are freely available on the Internet. Furthermore, they’ve encoded some of the logic right into the data (Math.Cos) and the implementation is incomplete, as it leaves us with a bunch of terms that we still have to do further processing on (albeit only a simple addition). Finally, we end up with different methods for each of the 6 variants of VSOP87, so we need something like 6 x 3 x 6 = 108 distinct methods for a full implementation.

Looking at AA+ and AAPlus, these take a similar approach, but they haven’t encoded the Cosine computation right into the data. However, consider this code snippet from AAPlus:

   1: public class Earth
   2:     {
   3:         public Earth() {}
   5:         static readonly VSOP87Coefficient [] g_L0EarthCoefficients = new VSOP87Coefficient[]
   6:         {
   7:         new VSOP87Coefficient(175347046, 0,         0 ),
   8:         new VSOP87Coefficient(3341656,   4.6692568,  6283.0758500 ),
   9:         new VSOP87Coefficient(34894,     4.62610,   12566.15170 ),
  10:         new VSOP87Coefficient(3497,      2.7441,    5753.3849 ),
  11:         '...3,000 or so lines elided for clarity
  12:         };
  14:         public static double EclipticLongitude(double JD)
  15:         {
  16:             double rho = (JD - 2451545) / 365250;
  17:             double rhosquared = rho*rho;
  18:             double rhocubed = rhosquared*rho;
  19:             double rho4 = rhocubed*rho;
  20:             double rho5 = rho4*rho;
  22:             //Calculate L0
  23:             int nL0Coefficients = g_L0EarthCoefficients.Length;
  24:             double L0 = 0;
  25:             for (int i=0; i<nL0Coefficients; i++)
  26:                 L0 += g_L0EarthCoefficients[i].A * Math.Cos(g_L0EarthCoefficients[i].B + g_L0EarthCoefficients[i].C*rho);
  28:             //Calculate L1
  29:             int nL1Coefficients = g_L1EarthCoefficients.Length;
  30:             double L1 = 0;
  31:             for (int i=0; i<nL1Coefficients; i++)
  32:                 L1 += g_L1EarthCoefficients[i].A * Math.Cos(g_L1EarthCoefficients[i].B + g_L1EarthCoefficients[i].C*rho);
  34:             //Calculate L2
  35:             int nL2Coefficients = g_L2EarthCoefficients.Length;
  36:             double L2 = 0;
  37:             for (int i=0; i<nL2Coefficients; i++)
  38:                 L2 += g_L2EarthCoefficients[i].A * Math.Cos(g_L2EarthCoefficients[i].B + g_L2EarthCoefficients[i].C*rho);
  40:             //Calculate L3
  41:             int nL3Coefficients = g_L3EarthCoefficients.Length;
  42:             double L3 = 0;
  43:             for (int i=0; i<nL3Coefficients; i++)
  44:                 L3 += g_L3EarthCoefficients[i].A * Math.Cos(g_L3EarthCoefficients[i].B + g_L3EarthCoefficients[i].C*rho);
  46:             //Calculate L4
  47:             int nL4Coefficients = g_L4EarthCoefficients.Length;
  48:             double L4 = 0;
  49:             for (int i=0; i<nL4Coefficients; i++)
  50:                 L4 += g_L4EarthCoefficients[i].A * Math.Cos(g_L4EarthCoefficients[i].B + g_L4EarthCoefficients[i].C*rho);
  52:             //Calculate L5
  53:             int nL5Coefficients = g_L5EarthCoefficients.Length;
  54:             double L5 = 0;
  55:             for (int i=0; i<nL5Coefficients; i++)
  56:                 L5 += g_L5EarthCoefficients[i].A * Math.Cos(g_L5EarthCoefficients[i].B + g_L5EarthCoefficients[i].C*rho);
  58:             double val = (L0 + L1*rho + L2*rhosquared + L3*rhocubed + L4*rho4 + L5*rho5) / 100000000.0;
  60:             //convert results back to degrees
  61:             val = CoordinateTransformation.MapTo0To360Range(CoordinateTransformation.RadiansToDegrees(val));
  62:             return val;
  63:         }


OK, so this code is computing the ecliptic longitude of the Earth. For each of the 6 terms (L0 to L5 in the code) the same computation is repeated over and over again. It is repeated another 6 times in each of the following methods:

   1: public static double EclipticLatitude(double JD) {/*...*/}
   2: public static double RadiusVector(double JD) {/*...*/}
   3: public static double EclipticLongitudeJ2000(double JD) {/*...*/}
   4: public static double EclipticLatitudeJ2000(double JD) {/*...*/}

If your eyes are not already bleeding from DRY violation, consider that all the above repetition is raised to another power, because each planet has its own class, each with its own copy of the VSOP87 code with boilerplate duplicate code in each class. This, and other libraries are full of SOLID violations such as this. They are full of static ‘helper methods’ and just not well structured as an object oriented design. Of course, the have classes and properties and methods, but simply having those things doesn’t make the code object oriented. The concepts need to be teased out and put back together using more modern approaches. Most of these libraries were ported from older languages and they just copied lines of code, with slight syntax tweaks to make it compile, without any real attempt at OO design.

VSOP87, as with a lot of astrophysics, is quite a complicated thing to understand, so the authors of these libraries can be forgiven to a large extent – I have spent many hours trying to really understand VSOP87 and although I still don’t fully grok the theory, I think I have a pretty good handle on how to implement the algorithm. There is a *lot* of improvement that can be done on all this code. I’m positive we can get much DRYer, we can improve the cohesion and eliminate all the close coupling and make things testable. We can split the VSOP87 code out into an ‘orbit engine’ and make it pluggable so we can plug in different orbit engines; we can pull the data back out of the code into data files or at least make it so we can somehow plug in new data if we need to.

Making a Start

OK, time to put my keyboard where my mouth is! I’ve written my first specification. It is a humble beginning, but it is a beginning. I was unsure where to start so I began with something I knew I could verify. I used the VSOP87 code generator to generate a set of methods for VSOP87B, and I picked just the L0 term. In my production code, I cheated by returning a constant (in TDD, all sins are forgiven to get to a green bar). I have not yet started to refactor the code with a proper implementation, that will be the subject of my next article. Here is the code so far:

   1: using Machine.Specifications;
   2: using TA.Orbits.ReferenceData;
   4: namespace TA.Orbits.Specifications
   5:     {
   6:     /*
   7:      * Calculate the position in space of the Earth relative to the Sun for a given date, time
   8:      * Give the answer in both cartesian coordinates (X,Y,Z)
   9:      * and sperical coordinates (Latitude, Longitude and Radius).
  10:      * 
  11:      * Use a reference implementation to verify the results.
  12:      */
  14:     [Subject(typeof(VSOP87B_EarthPositionSpherical), "reference data")]
  15:     public class when_computing_vsop87_term_l0_for_earth
  16:         {
  17:         Establish context = () =>
  18:             {
  19:             var rho = (J2000 - 2451545.0)/365250.0;
  20:             ReferenceEarthL0 = VSOP87B_EarthPositionSpherical.Earth_L0(rho);
  21:             OrbitEngine = new OrbitEngine.Vsop87.OrbitEngine();
  22:             };
  24:         Because of = () => ComputedEarthL0 = OrbitEngine.ComputeL0(J2000);
  26:         It should_compute_l0_to_match_the_reference_implementation =
  27:             () => ComputedEarthL0.ShouldBeCloseTo(ReferenceEarthL0, Tolerance);
  29:         static double ReferenceEarthL0;
  30:         static double ComputedEarthL0;
  31:         static OrbitEngine.Vsop87.OrbitEngine OrbitEngine;
  32:         static double Tolerance = 0.00000000000005; // 10E-14
  33:         const double J2000 = 2451545.0; // Julian date JD 2000.0
  34:         }
  35:     }


   1: namespace TA.OrbitEngine.Vsop87
   2:     {
   3:     public class OrbitEngine
   4:         {
   5:         public double ComputeL0(double julianDate)
   6:             {
   7:             return 1.75192386367183;
   8:             }
   9:         }
  10:     }


Humble beginnings, to be sure. I’m keeping my source code in a Git repository at http://stash.teamserver.tigranetworks.co.uk/projects/TNO/repos/ta.orbits/browse. In the next episode I’ll start to develop my code. Meanwhile, having stuck my head above the parapet by criticizing other astronomy libraries, I realize that I have opened myself to criticism in return. I welcome that. Please share your thoughts with me. You may well have better ideas than me and as long as you’re constructive, I’ll listen and learn from you.

1 comment(s)