Object Oriented Astronomy Part 2

Monday February 03 2014 - , , , , , , ,

Previously: http://www.tigranetworks.co.uk/blogs/electricdreams/object-oriented-astronomy-ndash-part-1/

Prologue

So here we are in the second installment, time to write some code! I’ve given myself the mission, for the moment, to compute the position of Earth in space (relative to the Sun) for a given date and time and to get the result in both rectangular (x, y, z) and spherical (latitude, longitude and radius) coordinates. I’m using an algorithm known as VSOP87 to do that. I want to be able to run everything from the original data files, so they can be updated in the future without changing any code and to check everything against a reference implementation. At the end of the last installment, I had written a single unit test ‘when_computing_vsop87_term_l0_for_earth’ and the implementation simply returned a fudged constant, getting me quickly to a green bar. I had already imported some reference data from the VSOP87 code generator web site. I now continue with the ‘refactor’ part of that test.

Refactor

First, a class to represent each individual term of a VSOP87 series:

   1: public sealed class Vsop87Term
   2:     {
   3:     /// <summary>
   4:     ///   Initializes a new instance of the <see cref="Vsop87Term" /> class.
   5:     ///   This constructor is appropriate when only the Amplitude A, Phase B and Frequency C
   6:     ///   values are required.
   7:     /// </summary>
   8:     /// <param name="amplitude">The amplitude.</param>
   9:     /// <param name="phase">The phase.</param>
  10:     /// <param name="frequency">The frequency.</param>
  11:     public Vsop87Term(double amplitude, double phase, double frequency)
  12:         {
  13:         AmplitudeA = amplitude;
  14:         PhaseB = phase;
  15:         FrequencyC = frequency;
  16:         }
  17:  
  18:     /// <summary>
  19:     ///   Gets the Amplitude (A) coefficient.
  20:     /// </summary>
  21:     /// <value>
  22:     ///   The Amplitude (A) coefficient expressed in au/(tjy**alpha) for the distances
  23:     ///   and in rad/(tjy**alpha) for angular variables.
  24:     /// </value>
  25:     public double AmplitudeA { get; private set; }
  26:  
  27:     /// <summary>
  28:     ///   Gets the Phase (B) coefficient.
  29:     /// </summary>
  30:     /// <value>The Phase (B) coefficient expressed in radians.</value>
  31:     public double PhaseB { get; private set; }
  32:  
  33:     /// <summary>
  34:     ///   Gets the Frequency (C) coefficient.
  35:     /// </summary>
  36:     /// <value>The Frequency (C) coefficient expressed in rad/tjy.</value>
  37:     public double FrequencyC { get; private set; }
  38:     };

This is a bit of a simplification as each row of data has many more terms than this, but we only need the last 3 terms of each row for spherical and rectangular coordinates.

Now, by hand-editing the VSOP87 data file, with plenty of use of column edit, it is a fairly simple matter to produce hard coded data, the beginning of which looks like this…

   1: using System.Collections.Generic;
   2: using TA.OrbitEngine.Vsop87;
   3:  
   4: internal static class Vsop87Data
   5:     {
   6:     internal static IList<Vsop87Term> Vsop87B_Earth_L0 = new List<Vsop87Term>
   7:         {
   8:         new Vsop87Term(1.75347045673, 0.00000000000, 0.00000000000),
   9:         new Vsop87Term(0.03341656453, 4.66925680415, 6283.07584999140),
  10:         new Vsop87Term(0.00034894275, 4.62610242189, 12566.15169998280),
  11:         // etc...

There are some 622 terms in that first series alone, with 6 such series for each coordinate variable and 3 coordinate variables. Although the length of the series does get dramatically shorter as the power increases, nevertheless there are well over 1,000 rows of data to be processed for each of the 3 coordinate variables. As the code generator web site states, entering that much data is indeed prone to error and for that reason, I plan to eventually use the original, unedited data files. For now though, hard coded data will do.

The formula to crunch all this data is relatively simple:

   1: public class OrbitEngine
   2:     {
   3:     public double ComputeL0(double julianDate)
   4:         {
   5:         //return 1.75192386367183;
   6:         // Iteratively apply the formula Tn = AT^alpha Cos(B + CT)
   7:         // Sum all Tn and return the sum.
   8:         var alpha = 0;  // Power of the term being computed
   9:         var thousandsOfJulianDays = (julianDate - 2451545.0) / 365250.0; // Thousands of Julian Days since JD2000.0
  10:         var tjdPowerAlpha = Math.Pow(thousandsOfJulianDays, alpha);
  11:         double sum = 0.0;
  12:         foreach (var term in Vsop87Data.Vsop87B_Earth_L0)
  13:             {
  14:             sum += term.AmplitudeA*tjdPowerAlpha*Math.Cos(term.PhaseB + term.FrequencyC*thousandsOfJulianDays);
  15:             }
  16:         return sum;
  17:         }
  18:     }

… and we have a green bar once again, so we know that the computation is good and matches the reference implementation supplied by the code generator web site. But we have only computed a single power term, for α=0 in the formula:

We need to compute the other terms for α=1, 2, 3, 4 and 5, so our computation needs to be generalized to accept the power (value of α) and also we need to be able to pass in the correct data. Our specification changes to:

   1: Because of = () => ComputedEarthL0 = OrbitEngine.ComputeVsop87Term(J2000, 0.0, Vsop87Data.Vsop87B_Earth_L0);

…and the computation method to:

   1: public double ComputeVsop87Term(double julianDate, double alpha, IList<Vsop87Term> series)
   2:     {
   3:     //return 1.75192386367183;
   4:     // Iteratively apply the formula Tn = AT^alpha Cos(B + CT)
   5:     // Sum all Tn and return the sum.
   6:     var thousandsOfJulianDays = (julianDate - 2451545.0)/365250.0; // Thousands of Julian Days since JD2000.0
   7:     var tjdPowerAlpha = Math.Pow(thousandsOfJulianDays, alpha);
   8:     var sum = 0.0;
   9:     foreach (var term in series)
  10:         sum += term.AmplitudeA*tjdPowerAlpha*Math.Cos(term.PhaseB + term.FrequencyC*thousandsOfJulianDays);
  11:     return sum;
  12:     }

Still green, we didn’t break anything. Now, it would be nice to check this with a different dataset, so I have imported the data for the L5 term (alpha=5) and added a new test for that. There is commonality between the two tests classes so I have factored out a context base class, as follows:

   1: using Machine.Specifications;
   2: using TA.Orbits.ReferenceData;
   3:  
   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:      */
  13:  
  14:     [Subject(typeof(VSOP87B_EarthPositionSpherical), "reference data")]
  15:     public class when_computing_vsop87_term_l0_for_earth : with_target_date_2014_jan_29_midday
  16:         {
  17:         Establish context = () => ReferenceTerm = VSOP87B_EarthPositionSpherical.Earth_L0(Rho);
  18:         Because of = () => ComputedTerm = OrbitEngine.ComputeVsop87Term(TargetDate, 0.0, Vsop87Data.Vsop87B_Earth_L0);
  19:  
  20:         It should_match_the_reference_implementation =
  21:             () => ComputedTerm.ShouldBeCloseTo(ReferenceTerm, Tolerance);
  22:         }
  23:  
  24:     [Subject(typeof(VSOP87B_EarthPositionSpherical), "reference data")]
  25:     public class when_computing_vsop87_term_l5_for_earth : with_target_date_2014_jan_29_midday
  26:         {
  27:         Establish context = () => ReferenceTerm = VSOP87B_EarthPositionSpherical.Earth_L5(Rho);
  28:         Because of = () => ComputedTerm = OrbitEngine.ComputeVsop87Term(TargetDate, 5.0, Vsop87Data.Vsop87B_Earth_L5);
  29:         It should_match_the_reference_implementation =
  30:             () => ComputedTerm.ShouldBeCloseTo(ReferenceTerm, Tolerance);
  31:         }
  32:  
  33:     public class with_target_date_2014_jan_29_midday
  34:         {
  35:         protected const double J2000 = 2451545.0; // Julian date JD 2000.0
  36:         protected const double Tolerance = 0.00000000000005; // 10E-14
  37:         protected static OrbitEngine.Vsop87.OrbitEngine OrbitEngine;
  38:         protected static double Rho;
  39:         protected static double TargetDate;
  40:         protected static double ReferenceTerm;
  41:         protected static double ComputedTerm;
  42:         Establish context = () =>
  43:             {
  44:             TargetDate = 2456687; // 29/Jan/2014 midday
  45:             Rho = (TargetDate - 2451545.0)/365250.0;
  46:             OrbitEngine = new OrbitEngine.Vsop87.OrbitEngine();
  47:             };
  48:         }
  49:     }

The data for the L5 term is short, here it is in its entirety:

   1: internal static IList<Vsop87Term> Vsop87B_Earth_L5 = new List<Vsop87Term>
   2:     {
   3:     // VSOP87 VERSION B2    EARTH     VARIABLE 1 (LBR)       *T**5      4 TERMS    HELIOCENTRIC DYNAMICAL ECLIPTIC AND EQUINOX J2000      
   4:     new Vsop87Term(0.00000000172, 2.74854172392, 6283.07584999140),
   5:     new Vsop87Term(0.00000000050, 2.01352986713, 155.42039943420),
   6:     new Vsop87Term(0.00000000028, 2.93369985477, 12566.15169998280),
   7:     new Vsop87Term(0.00000000005, 1.93829214518, 18849.22754997420),
   8:     };

The reference data we imported from the code generator web site already contains this data, so our new test runs and passes without any issues.

image

This is the state of the code at http://stash.teamserver.tigranetworks.co.uk/projects/TNO/repos/ta.orbits/commits/08967f6026cbb7d1aaf747dcbdf3b909379e423b

Moving on, we’ve only calculated two of the six terms required for the latitude coordinate, the next step seems obvious: compute all siz terms and return the latitude. Here’s the spec:

   1: [Subject(typeof(VSOP87B_EarthPositionSpherical), "reference data")]
   2: public class when_computing_vsop87_latitude_for_earth : with_target_date_2014_jan_29_midday
   3:     {
   4:     Establish context = () =>
   5:             ReferenceLatitude =
   6:                 VSOP87B_EarthPositionSpherical.Earth_L0(Rho) + VSOP87B_EarthPositionSpherical.Earth_L1(Rho) +
   7:                 VSOP87B_EarthPositionSpherical.Earth_L2(Rho) + VSOP87B_EarthPositionSpherical.Earth_L3(Rho) +
   8:                 VSOP87B_EarthPositionSpherical.Earth_L4(Rho) + VSOP87B_EarthPositionSpherical.Earth_L5(Rho);
   9:     Because of = () => ComputedLatitude = OrbitEngine.ComputeVsop87Series(TargetDate, Vsop87Data.Vsop87B_Earth_Latitude);
  10:     It should_match_the_reference_implementation =
  11:         () => ComputedLatitude.ShouldBeCloseTo(ReferenceLatitude, Tolerance);
  12: }

…and the method to perform the computation:

   1: public double ComputeVsop87Series(double targetDate, List<List<Vsop87Term>> seriesData)
   2:     {
   3:     var alpha = 0;  // series power
   4:     var sum = 0.0;
   5:     foreach (var term in seriesData)
   6:         {
   7:         sum += ComputeVsop87Term(targetDate, alpha, term);
   8:         ++alpha;
   9:         }
  10:     return sum;
  11:     }

Simple enough, and it passes, meaning that we still agree with the reference implementation.

Fun and Games With Data Files

At this point, I could continue hard coding the data for the longitude and radius variables, but the objective is to be able to work from the original, unedited data files, so the next step seems to be to load the data from a file. I specified it like this:

   1: [Subject(typeof(VSOP87B_EarthPositionSpherical), "reference data")]
   2: public class when_reading_earth_vsop87_data_from_a_text_file
   3: {
   4:     protected const double Tolerance = 0.00000000001;
   5:     Because of = () => EarthData = Vsop87DataReader.LoadVsop87DataFromFile("VSOP87B.ear");
   6:     It should_be_the_same_as_the_hard_coded_data = () => EarthData.VariableData['L'].ShouldBeLikeObjectGraph(Vsop87Data.Vsop87B_Earth_Latitude, Tolerance) ;
   7:     static Vsop87Solution EarthData;
   8: }

and the stubbed out implementation, giving a failing test:

   1: public object LoadVsop87DataFromFile(string filename)
   2:     {
   3:     throw new NotImplementedException();
   4:     }

Now, I’m writing this blog article retrospectively and this test led to some fun and games. I felt like smaller steps were needed here, but I couldn’t work out how to achieve that other than to go for what I shall optimistically refer to as an ‘obvious implementation’. Ahem. I actually spent several hours with a red test bar, which felt completely wrong. The code I came up with probably reflects this and I’d be interested in hearing how others would have dealt with it, or indeed if you can come up with a better implementation.

My first step was to come up with a couple of regular expressions for recognizing a header row and a data row withing the VSOP87 data files. I used Roy Osherove’s Regulator to produce the expressions using actual test data from the VSOP87 data files. Here’s what I came up with:

   1: const string Vsop87HeaderRegexPattern =
   2:     @"^\s*VSOP87 VERSION (?<VsopVersion>\w+)\s+(?<Body>\w+)\s+VARIABLE (?<Variable>\d+)\s+\((?<Variables>\w+)\)\s+\*T\*\*(?<Power>\d+)\s+(?<Terms>\d+)\s+TERMS\s+(?<Description>.*)";
   3: const string Vsop87DataRowRegexPattern = @"^(\s+[0-9.]+){16}\s+(?<A>[0-9.]+)\s+(?<B>[0-9.]+)\s+(?<C>[0-9.]+).*$";

Then I added the whole of the file reading implementation in one go, along with a class Vsop87Solution to hold the results. Here is the whole lump of code:

   1: public class Vsop87DataReader
   2:     {
   3:     const string Vsop87HeaderRegexPattern =
   4:         @"^\s*VSOP87 VERSION (?<VsopVersion>\w+)\s+(?<Body>\w+)\s+VARIABLE (?<Variable>\d+)\s+\((?<Variables>\w+)\)\s+\*T\*\*(?<Power>\d+)\s+(?<Terms>\d+)\s+TERMS\s+(?<Description>.*)";
   5:     const string Vsop87DataRowRegexPattern = @"^(\s*[+-]*[0-9.]+){16}\s+(?<A>[0-9.]+)\s+(?<B>[0-9.]+)\s+(?<C>[0-9.]+).*$";
   6:     internal static Regex Vsop87HeaderRegex = new Regex(Vsop87HeaderRegexPattern, RegexOptions.Compiled);
   7:     internal static Regex Vsop87DataRowRegex = new Regex(Vsop87DataRowRegexPattern, RegexOptions.Compiled);
   8:  
   9:     public static Vsop87Solution LoadVsop87DataFromFile(string filename)
  10:         {
  11:         Diagnostics.TraceInfo("Loading VSOP87 solution from {0}", filename);
  12:         var dataFile = Path.Combine(@".\Vsop87_data", filename);
  13:         var inputStream = File.OpenText(dataFile);
  14:         var header = inputStream.ReadLine();
  15:         var currentLine = header;
  16:         var solution = Vsop87Solution.FromHeaderString(header);
  17:         foreach (char variable in solution.Variables)
  18:             {
  19:             Diagnostics.TraceVerbose("Loading variable {0}", variable);
  20:             var variableIndex = solution.Variables.IndexOf(variable);
  21:             var series = LoadVariableSeries(solution, variableIndex, inputStream);
  22:             solution.VariableData[variable] = series;
  23:             }
  24:         return solution;
  25:         }
  26:  
  27:     /// <summary>
  28:     /// Loads the variable series.
  29:     /// </summary>
  30:     /// <param name="solution">The solution.</param>
  31:     /// <param name="variableIndex">Index of the variable.</param>
  32:     /// <returns>IEnumerable{IEnumerable{Vsop87Term}}.</returns>
  33:     static IEnumerable<IEnumerable<Vsop87Term>> LoadVariableSeries(Vsop87Solution solution, int variableIndex, StreamReader inputStream)
  34:         {
  35:         var majorSeries = new List<List<Vsop87Term>>();
  36:         int power = 0;
  37:         string currentVariable = (++variableIndex).ToString(CultureInfo.InvariantCulture);
  38:         string headerVariable = currentVariable;
  39:         do
  40:             {
  41:             var minorSeries = new List<Vsop87Term>();
  42:             var header = LoadPowerSeries(minorSeries, inputStream);
  43:             majorSeries.Add(minorSeries);
  44:             ++power;
  45:  
  46:             if (string.IsNullOrEmpty(header))
  47:                 break;  // End of stream
  48:  
  49:             var match = Vsop87HeaderRegex.Match(header);
  50:             headerVariable = match.Groups["Variable"].Value;
  51:             }
  52:         while (headerVariable == currentVariable);
  53:         return majorSeries;
  54:         }
  55:  
  56:     /// <summary>
  57:     /// Loads the terms in a power series for one power of one coordinate variable,
  58:     /// up to the next header row or end-of-stream. If a header row is found,
  59:     /// it is returned.
  60:     /// </summary>
  61:     /// <param name="minorSeries">The minor series.</param>
  62:     /// <param name="inputStream">The input stream.</param>
  63:     /// <returns>System.String containing the next header, or String.Empty if at end-of-stream.</returns>
  64:     /// <exception cref="System.InvalidOperationException">Data file was not in the expected format (could not recognize a data row or a header row)</exception>
  65:     static string LoadPowerSeries(List<Vsop87Term> minorSeries, StreamReader inputStream)
  66:         {
  67:         while (!inputStream.EndOfStream)
  68:             {
  69:             var line = inputStream.ReadLine();
  70:             Diagnostics.TraceVerbose("Read line: {0}", line);
  71:             var rowMatch = Vsop87DataRowRegex.Match(line);
  72:             if (rowMatch.Success)
  73:                 {
  74:                 var termA = double.Parse(rowMatch.Groups["A"].Value, CultureInfo.InvariantCulture);
  75:                 var termB = double.Parse(rowMatch.Groups["B"].Value, CultureInfo.InvariantCulture);
  76:                 var termC = double.Parse(rowMatch.Groups["C"].Value, CultureInfo.InvariantCulture);
  77:                 var term = new Vsop87Term(termA, termB, termC);
  78:                 Diagnostics.TraceVerbose("Recognized data row, A={0} B={1} C={2}", termA, termB, termC);
  79:                 minorSeries.Add(term);
  80:                 continue;
  81:                 }
  82:             var headerMatch = Vsop87HeaderRegex.Match(line);
  83:             if (headerMatch.Success)
  84:                 {
  85:                 Diagnostics.TraceVerbose("Recognized header row");
  86:                 return line;
  87:                 }
  88:             throw new InvalidOperationException("Data file was not in the expected format (could not recognize a data row or a header row)");
  89:             }
  90:         Diagnostics.TraceVerbose("Encountered end-of-stream while reading power series terms");
  91:         return string.Empty;    // At end of stream, no header.
  92:         }
  93:     }
  94:  
  95: /// <summary>
  96: /// Class Vsop87Solution. Holds all of the data for one body and one VSOP87 version
  97: /// (coordinate system and reference frame). VSOP87 has 6 distinct variants so there
  98: /// may be up to 6 separate solutions for each body.
  99: /// </summary>
 100: public class Vsop87Solution
 101:     {
 102:         /// <summary>
 103:         /// Gets the name of the body to which this solution applies.
 104:         /// </summary>
 105:         /// <value>The body name as an upper-case string.</value>
 106:     public string Body { get; private set; }
 107:     /// <summary>
 108:     /// Gets the VSOP87 version (variant) represented by this data.
 109:     /// A VSOP87 version defines the coordinate system and reference frame,
 110:     /// e.g. HELIOCENTRIC DYNAMICAL ECLIPTIC and EQUINOX J2000
 111:     /// </summary>
 112:     /// <value>The version, a string containing a single upper-case letter and a single decimal digit.</value>
 113:     public string Version { get; private set; }
 114:     /// <summary>
 115:     /// Gets the coordinate variables that can be computed with this data.
 116:     /// </summary>
 117:     /// <value>The variables, one character for each variable, concatenated together in a string.</value>
 118:     public string Variables { get; private set; }
 119:     /// <summary>
 120:     /// Gets the description of the coordinate system and reference frame represented by this data.
 121:     /// </summary>
 122:     /// <value>The description, as obtained from the VSOP87 data file.</value>
 123:     public string Description { get; private set; }
 124:  
 125:     /// <summary>
 126:     /// Gets the variable data for this VSOP87 solution.
 127:     /// The variables (coordinates) are contained in a dictionary, where the key
 128:     /// is the variable name and the value is the series data for that variable.
 129:     /// </summary>
 130:     /// <value>The variable data.</value>
 131:     public IDictionary<char, IEnumerable<IEnumerable<Vsop87Term>>> VariableData { get; private set; }
 132:  
 133:     public Vsop87Solution()
 134:         {
 135:         VariableData = new Dictionary<char, IEnumerable<IEnumerable<Vsop87Term>>>();
 136:         }
 137:     internal static Vsop87Solution FromHeaderString(string header)
 138:         {
 139:         var matches = Vsop87DataReader.Vsop87HeaderRegex.Match(header);
 140:         if (!matches.Success)
 141:             throw new ArgumentException("The supplied header was invalid", "header");
 142:         var body = matches.Groups["Body"].Value;
 143:         var version = matches.Groups["VsopVersion"].Value;
 144:         var variables = matches.Groups["Variables"].Value;
 145:         var description = matches.Groups["Description"].Value;
 146:         Diagnostics.TraceInfo("Creating VSOP87 solution with Body={0} Version={1} Variables={2} Description={3}",
 147:             body,
 148:             version,
 149:             variables,
 150:             description);
 151:        return new Vsop87Solution {Body = body, Version = version, Variables = variables, Description = description};
 152:         }
 153:     }

The eagle-eyed amongst you will have spotted the strange ShouldBeLikeObjectGraph() assertion in my test. I originally tried to use ShouldBeLike() to compare the two object graphs. This *almost* succeeded except that, out of the 1,000+ rows of data in the latitude variable, the comparison failed on about 4 of the entries. How annoying! This was down to the imprecise nature of floating point arithmetic and ShouldBeLike() uses an equality comparison regardless of the data type, with no facility to specify a tolerance as is the case with ShouldBeCloseTo(). Interestingly, for some values, double.Parse() produces a different value than an object initializer expression given the same literal text; this screen snip taken during debugging the issue illustrates the problem nicely – the difference is almost imperceptibly small yet a difference it is, and the test fails.

image

I worked around this by making my own version of ShouldBeLike() called ShouldBeLikeObjectGraph() that accepts a Tolerance parameter and uses it when comparing doubles.

image

It worked, but it feels like a bit of a hack. Later on, this test became redundant along with all the hard coded data and I was able to remove it, but at this point I needed to prove that my data file loading class produced the same object graph as the hard coded data. I will not reproduce the code for the ShouldBeLikeObjectGraph() method here but if you are interested you can find it in the repository in this commit:
http://stash.teamserver.tigranetworks.co.uk/projects/TNO/repos/ta.orbits/commits/664ccfd3167cb0c29f8afc65c55d3ed762492787 (note that this test later becomes redundant as is removed).

This led on to the first test where all of the data is loaded from a data file and returns a fully computed coordinate variable, radius in this case:

   1: [Subject(typeof(Vsop87OrbitEngine), "data loaded from file")]
   2: public class when_computing_vsop87_radius_variable_for_earth_with_data_loaded_from_file : with_target_date_2014_jan_29_midday
   3:     {
   4:     Establish context = () =>
   5:         {
   6:         ReferenceRadius =
   7:             VSOP87B_EarthPositionSpherical.Earth_R0(Rho) + VSOP87B_EarthPositionSpherical.Earth_R1(Rho) +
   8:             VSOP87B_EarthPositionSpherical.Earth_R2(Rho) + VSOP87B_EarthPositionSpherical.Earth_R3(Rho) +
   9:             VSOP87B_EarthPositionSpherical.Earth_R4(Rho) + VSOP87B_EarthPositionSpherical.Earth_R5(Rho);
  10:         EarthData = Vsop87DataReader.LoadVsop87DataFromFile("VSOP87B.ear");
  11:         };
  12:     Because of = () => ComputedRadius = Vsop87OrbitEngine.ComputeVsop87Series(TargetDate, EarthData.VariableData['R']);
  13:     It should_match_the_reference_implementation =
  14:         () => ComputedRadius.ShouldBeCloseTo(ReferenceRadius, Tolerance);
  15:     static Vsop87Solution EarthData;
  16:     }

and the results:

image

In this last test, the name of the data file is specified directly. I didn’t like exposing this detail to the user, so I decided to introduce a method that would automatically select the correct data file. The user specifies the solar system body, coordinate system and reference frame required and the method works out which VSOP87 variant is needed and builds the correct file name.

   1: [Subject(typeof(Vsop87DataReader), "Data file selection")]
   2: public class when_selecting_a_data_file
   3:     {
   4:     It should_be_variant_base_emb_for_earth_moon_barycentre_j2000_elliptic_elements = () =>
   5:         Vsop87DataReader.SelectDataFile(SolarSystemBody.EarthMoonBarycentre,
   6:             CoordinateSystem.HeliocentricEllipticElements,
   7:             ReferenceFrame.EquinoxJ2000).ShouldEqual("VSOP87.emb");
   8:     It should_be_variant_a_mer_for_mercury_j2000_rectangular = () =>
   9:         Vsop87DataReader.SelectDataFile(SolarSystemBody.Mercury,
  10:             CoordinateSystem.HeliocentricRectangularCoordinates,
  11:             ReferenceFrame.EquinoxJ2000).ShouldEqual("VSOP87A.mer");
  12:     It should_be_variant_b_ven_for_venus_j2000_spherical = () =>
  13:         Vsop87DataReader.SelectDataFile(SolarSystemBody.Venus,
  14:             CoordinateSystem.HeliocentricSphericalCoordinates,
  15:             ReferenceFrame.EquinoxJ2000).ShouldEqual("VSOP87B.ven");
  16:     It should_be_variant_c_nep_for_neptune_jnow_rectangular = () =>
  17:         Vsop87DataReader.SelectDataFile(SolarSystemBody.Neptune,
  18:             CoordinateSystem.HeliocentricRectangularCoordinates,
  19:             ReferenceFrame.EquinoxJNow).ShouldEqual("VSOP87C.nep");
  20:     It should_be_variant_d_ear_for_earth_jnow_spherical = () =>
  21:         Vsop87DataReader.SelectDataFile(SolarSystemBody.Earth,
  22:             CoordinateSystem.HeliocentricSphericalCoordinates,
  23:             ReferenceFrame.EquinoxJNow).ShouldEqual("VSOP87D.ear");
  24:     It should_be_variant_e_sat_for_saturn_j2000_barycentric = () =>
  25:         Vsop87DataReader.SelectDataFile(SolarSystemBody.Saturn,
  26:             CoordinateSystem.BarycentricRectangularCoordinates,
  27:             ReferenceFrame.EquinoxJ2000).ShouldEqual("VSOP87E.sat");
  28:     }

This last test deviates from the normal Arrange/Act/Assert pattern but it seems to read well and make sense in that format. While implementing the file selection method, I realized that not all combinations were valid (i.e. VSOP87 did not include the data) and therefore I produced a series of specs to explicitly check for these invalid combinations:

   1: [Subject(typeof(Vsop87DataReader), "data file selection")]
   2: public class when_selecting_an_unsupported_configuration_of_barycentric_coordinates_and_jnow
   3: {
   4:     Because of = () => Thrown = Catch.Exception(() =>
   5:         Vsop87DataReader.SelectDataFile(SolarSystemBody.Mars,
   6:             CoordinateSystem.BarycentricRectangularCoordinates,
   7:             ReferenceFrame.EquinoxJNow));
   8:     It should_throw_not_supported_exception = () =>
   9:         Thrown.ShouldBeOfType<NotSupportedException>();
  10:     It should_contain_the_body =
  11:         () => Thrown.Data.ShouldContain(new DictionaryEntry("SolarSystemBody", SolarSystemBody.Mars));
  12:     It should_contain_the_coordinate_system =
  13:         () => Thrown.Data.ShouldContain(new DictionaryEntry("CoordinateSystem", CoordinateSystem.BarycentricRectangularCoordinates));
  14:     It should_contain_the_reference_frame =
  15:         () => Thrown.Data.ShouldContain(new DictionaryEntry("ReferenceFrame", ReferenceFrame.EquinoxJNow));
  16:     static Exception Thrown;
  17: }
  18:  
  19: [Subject(typeof(Vsop87DataReader), "data file selection")]
  20: public class when_selecting_an_unsupported_configuration_of_elliptic_elements_and_jnow
  21: {
  22:     Because of = () => Thrown = Catch.Exception(() =>
  23:         Vsop87DataReader.SelectDataFile(SolarSystemBody.Mars,
  24:             CoordinateSystem.HeliocentricEllipticElements,
  25:             ReferenceFrame.EquinoxJNow));
  26:     It should_throw_not_supported_exception = () =>
  27:         Thrown.ShouldBeOfType<NotSupportedException>();
  28:     It should_contain_the_body =
  29:         () => Thrown.Data.ShouldContain(new DictionaryEntry("SolarSystemBody", SolarSystemBody.Mars));
  30:     It should_contain_the_coordinate_system =
  31:         () => Thrown.Data.ShouldContain(new DictionaryEntry("CoordinateSystem", CoordinateSystem.HeliocentricEllipticElements));
  32:     It should_contain_the_reference_frame =
  33:         () => Thrown.Data.ShouldContain(new DictionaryEntry("ReferenceFrame", ReferenceFrame.EquinoxJNow));
  34:     static Exception Thrown;
  35: }
  36:  
  37: [Subject(typeof(Vsop87DataReader), "data file selection")]
  38: public class when_selecting_an_unsupported_configuration_of_the_sun_and_non_barycentric_coordinates
  39: {
  40:     Because of = () => Thrown = Catch.Exception(() =>
  41:         Vsop87DataReader.SelectDataFile(SolarSystemBody.Sun,
  42:             CoordinateSystem.HeliocentricRectangularCoordinates,
  43:             ReferenceFrame.EquinoxJ2000));
  44:     It should_throw_not_supported_exception = () =>
  45:         Thrown.ShouldBeOfType<NotSupportedException>();
  46:     It should_contain_the_body =
  47:         () => Thrown.Data.ShouldContain(new DictionaryEntry("SolarSystemBody", SolarSystemBody.Sun));
  48:     It should_contain_the_coordinate_system =
  49:         () => Thrown.Data.ShouldContain(new DictionaryEntry("CoordinateSystem", CoordinateSystem.HeliocentricRectangularCoordinates));
  50:     It should_contain_the_reference_frame =
  51:         () => Thrown.Data.ShouldContain(new DictionaryEntry("ReferenceFrame", ReferenceFrame.EquinoxJ2000));
  52:     static Exception Thrown;
  53: }
  54:  
  55: [Subject(typeof(Vsop87DataReader), "data file selection")]
  56: public class when_selecting_an_unsupported_configuration_of_earth_moon_barycentre_and_jnow
  57: {
  58:     Because of = () => Thrown = Catch.Exception(() =>
  59:         Vsop87DataReader.SelectDataFile(SolarSystemBody.EarthMoonBarycentre,
  60:             CoordinateSystem.HeliocentricRectangularCoordinates,
  61:             ReferenceFrame.EquinoxJNow));
  62:     It should_throw_not_supported_exception = () =>
  63:         Thrown.ShouldBeOfType<NotSupportedException>();
  64:     It should_contain_the_body =
  65:         () => Thrown.Data.ShouldContain(new DictionaryEntry("SolarSystemBody", SolarSystemBody.EarthMoonBarycentre));
  66:     It should_contain_the_coordinate_system =
  67:         () => Thrown.Data.ShouldContain(new DictionaryEntry("CoordinateSystem", CoordinateSystem.HeliocentricRectangularCoordinates));
  68:     It should_contain_the_reference_frame =
  69:         () => Thrown.Data.ShouldContain(new DictionaryEntry("ReferenceFrame", ReferenceFrame.EquinoxJNow));
  70:     static Exception Thrown;
  71: }
  72:  
  73: [Subject(typeof(Vsop87DataReader), "data file selection")]
  74: public class when_selecting_an_unsupported_configuration_of_earth_moon_barycentre_and_spherical_coordinates
  75: {
  76:     Because of = () => Thrown = Catch.Exception(() =>
  77:         Vsop87DataReader.SelectDataFile(SolarSystemBody.EarthMoonBarycentre,
  78:             CoordinateSystem.HeliocentricSphericalCoordinates,
  79:             ReferenceFrame.EquinoxJ2000));
  80:     It should_throw_not_supported_exception = () =>
  81:         Thrown.ShouldBeOfType<NotSupportedException>();
  82:     It should_contain_the_body =
  83:         () => Thrown.Data.ShouldContain(new DictionaryEntry("SolarSystemBody", SolarSystemBody.EarthMoonBarycentre));
  84:     It should_contain_the_coordinate_system =
  85:         () => Thrown.Data.ShouldContain(new DictionaryEntry("CoordinateSystem", CoordinateSystem.HeliocentricSphericalCoordinates));
  86:     It should_contain_the_reference_frame =
  87:         () => Thrown.Data.ShouldContain(new DictionaryEntry("ReferenceFrame", ReferenceFrame.EquinoxJ2000));
  88:     static Exception Thrown;
  89: }
  90:  
  91: [Subject(typeof(Vsop87DataReader), "data file selection")]
  92: public class when_selecting_an_unsupported_configuration_of_earth_moon_barycentre_and_barycentric_coordinates
  93: {
  94:     Because of = () => Thrown = Catch.Exception(() =>
  95:         Vsop87DataReader.SelectDataFile(SolarSystemBody.EarthMoonBarycentre,
  96:             CoordinateSystem.BarycentricRectangularCoordinates,
  97:             ReferenceFrame.EquinoxJ2000));
  98:     It should_throw_not_supported_exception = () =>
  99:         Thrown.ShouldBeOfType<NotSupportedException>();
 100:     It should_contain_the_body =
 101:         () => Thrown.Data.ShouldContain(new DictionaryEntry("SolarSystemBody", SolarSystemBody.EarthMoonBarycentre));
 102:     It should_contain_the_coordinate_system =
 103:         () => Thrown.Data.ShouldContain(new DictionaryEntry("CoordinateSystem", CoordinateSystem.BarycentricRectangularCoordinates));
 104:     It should_contain_the_reference_frame =
 105:         () => Thrown.Data.ShouldContain(new DictionaryEntry("ReferenceFrame", ReferenceFrame.EquinoxJ2000));
 106:     static Exception Thrown;
 107: }
 108:  
 109: [Subject(typeof(Vsop87DataReader), "data file selection")]
 110: public class when_selecting_an_unsupported_configuration_of_earth_and_elliptic_elements
 111: {
 112:     Because of = () => Thrown = Catch.Exception(() =>
 113:         Vsop87DataReader.SelectDataFile(SolarSystemBody.Earth,
 114:             CoordinateSystem.HeliocentricEllipticElements,
 115:             ReferenceFrame.EquinoxJ2000));
 116:     It should_throw_not_supported_exception = () =>
 117:         Thrown.ShouldBeOfType<NotSupportedException>();
 118:     It should_contain_the_body =
 119:         () => Thrown.Data.ShouldContain(new DictionaryEntry("SolarSystemBody", SolarSystemBody.Earth));
 120:     It should_contain_the_coordinate_system =
 121:         () => Thrown.Data.ShouldContain(new DictionaryEntry("CoordinateSystem", CoordinateSystem.HeliocentricEllipticElements));
 122:     It should_contain_the_reference_frame =
 123:         () => Thrown.Data.ShouldContain(new DictionaryEntry("ReferenceFrame", ReferenceFrame.EquinoxJ2000));
 124:     static Exception Thrown;
 125: }

…and here’s the method…

   1: /// <summary>
   2: /// Selects the VSOP87 data file based on the body required and the type of coordinates required.
   3: /// </summary>
   4: /// <param name="solarSystemBody">The solar system body.</param>
   5: /// <param name="coordinateSystem">The coordinate system.</param>
   6: /// <param name="referenceFrame">The reference frame.</param>
   7: /// <returns>A string containing the filename and extension (no path) of the VSOP87 data file to use.</returns>
   8: internal static string SelectDataFile(SolarSystemBody solarSystemBody, CoordinateSystem coordinateSystem, ReferenceFrame referenceFrame)
   9:     {
  10:     string body;
  11:     string variant;
  12:     const string fileNameFormat = "VSOP87{0}.{1}";
  13:     switch (solarSystemBody)
  14:         {
  15:         case SolarSystemBody.Sun:
  16:             // when the Sun is the target body only barycentric coordinates and equinox J2000 are valid.
  17:             if (coordinateSystem != CoordinateSystem.BarycentricRectangularCoordinates ||
  18:                 referenceFrame != ReferenceFrame.EquinoxJ2000)
  19:                 throw NotSupportedExceptionBuilder(
  20:                     "The Sun may only be used with barycentric coordinates and equinox J2000",
  21:                     solarSystemBody,
  22:                     coordinateSystem,
  23:                     referenceFrame);
  24:             body = "sun";
  25:             break;
  26:         case SolarSystemBody.Mercury:
  27:             body = "mer";
  28:             break;
  29:         case SolarSystemBody.Venus:
  30:             body = "ven";
  31:             break;
  32:         case SolarSystemBody.Earth:
  33:             // When Earth is the target body, elliptic elements are not valid.
  34:             if (coordinateSystem == CoordinateSystem.HeliocentricEllipticElements)
  35:                 throw NotSupportedExceptionBuilder(
  36:                     "Earth cannot be used with elliptic elements, please use rectangular or spherical coordinates instead.",
  37:                     solarSystemBody,
  38:                     coordinateSystem,
  39:                     referenceFrame);
  40:             body = "ear";
  41:             break;
  42:         case SolarSystemBody.Mars:
  43:             body = "mar";
  44:             break;
  45:         case SolarSystemBody.Jupiter:
  46:             body = "jup";
  47:             break;
  48:         case SolarSystemBody.Saturn:
  49:             body = "sat";
  50:             break;
  51:         case SolarSystemBody.Uranus:
  52:             body = "ura";
  53:             break;
  54:         case SolarSystemBody.Neptune:
  55:             body = "nep";
  56:             break;
  57:         case SolarSystemBody.EarthMoonBarycentre:
  58:             // When Earth-Moon Barycenter is the target body, only elliptical elements or rectangular J2000 coordinates are valid.
  59:             if (coordinateSystem == CoordinateSystem.HeliocentricSphericalCoordinates
  60:                 || referenceFrame == ReferenceFrame.EquinoxJNow ||
  61:                 coordinateSystem == CoordinateSystem.BarycentricRectangularCoordinates)
  62:                 throw NotSupportedExceptionBuilder(
  63:                     "Earth-Moon Barycentre is only valid with equinox J2000 and either elliptic or rectangular coordinates",
  64:                     solarSystemBody,
  65:                     coordinateSystem,
  66:                     referenceFrame);
  67:             body = "emb";
  68:             break;
  69:         default:
  70:             throw NotSupportedExceptionBuilder("The specified body is not supported by VSOP87",
  71:                 solarSystemBody,
  72:                 coordinateSystem,
  73:                 referenceFrame);
  74:         }
  75:  
  76:     switch (coordinateSystem)
  77:         {
  78:             case CoordinateSystem.HeliocentricEllipticElements:
  79:             if (referenceFrame == ReferenceFrame.EquinoxJ2000)
  80:                 {
  81:                 variant = string.Empty;
  82:                 }
  83:             else
  84:                 {
  85:                 throw NotSupportedExceptionBuilder(
  86:                     "Equinox of date (JNow) is not supported for elliptic elements, please use J2000, rectangular or spherical coordinates instead.",
  87:                     solarSystemBody,
  88:                     coordinateSystem,
  89:                     referenceFrame);
  90:                 }
  91:             break;
  92:         case CoordinateSystem.HeliocentricRectangularCoordinates:
  93:             variant = referenceFrame == ReferenceFrame.EquinoxJ2000 ? "A" : "C";
  94:             break;
  95:             case CoordinateSystem.HeliocentricSphericalCoordinates:
  96:             variant = referenceFrame == ReferenceFrame.EquinoxJ2000 ? "B" : "D";
  97:             break;
  98:             case CoordinateSystem.BarycentricRectangularCoordinates:
  99:             if (referenceFrame==ReferenceFrame.EquinoxJ2000)
 100:                 {
 101:                 variant = "E";
 102:                 }
 103:             else
 104:                 {
 105:                 throw NotSupportedExceptionBuilder(
 106:                     "Equinox of date (JNow) is not supported for barycentric coordinates, please use J2000 or heliocentric coordinates instead.",
 107:                     solarSystemBody,
 108:                     coordinateSystem,
 109:                     referenceFrame);
 110:                 }
 111:             break;
 112:         default:
 113:             throw NotSupportedExceptionBuilder("The slected coordinate system is not supported by VSOP87",
 114:                 solarSystemBody,
 115:                 coordinateSystem,
 116:                 referenceFrame);
 117:         }
 118:  
 119:     var filename = string.Format(fileNameFormat, variant, body);
 120:     return filename;
 121:     }
 122:  
 123: static Exception NotSupportedExceptionBuilder(string message, SolarSystemBody solarSystemBody, CoordinateSystem coordinateSystem, ReferenceFrame referenceFrame)
 124:     {
 125:     var ex = new NotSupportedException(message);
 126:     ex.Data.Add("SolarSystemBody", solarSystemBody);
 127:     ex.Data.Add("CoordinateSystem", coordinateSystem);
 128:     ex.Data.Add("ReferenceFrame", referenceFrame);
 129:     return ex;
 130:     }
 131: }

image

Bringing It All Together

We’re almost there now. All of the plumbing is in place. In the next phase, I add classes to represent rectangular and spherical coordinates which are the results we wanted. I implement a pair of methods to compute those coordinates, which now completely hide all of the details of the VSOP87 implementation from the user. Here are the tests and the methods:

   1: [Subject(typeof(Vsop87OrbitEngine), "compute coordinates")]
   2: public class when_computing_spherical_j2000_coordinates_for_earth : with_target_date_2014_jan_29_midday
   3: {
   4:     Establish context = () =>
   5:     {
   6:         var latitude = VSOP87B_EarthPositionSphericalJ2000.Earth_L0(Rho) + VSOP87B_EarthPositionSphericalJ2000.Earth_L1(Rho) +
   7:                        VSOP87B_EarthPositionSphericalJ2000.Earth_L2(Rho) + VSOP87B_EarthPositionSphericalJ2000.Earth_L3(Rho) +
   8:                        VSOP87B_EarthPositionSphericalJ2000.Earth_L4(Rho) + VSOP87B_EarthPositionSphericalJ2000.Earth_L5(Rho);
   9:         var longitude = VSOP87B_EarthPositionSphericalJ2000.Earth_B0(Rho) + VSOP87B_EarthPositionSphericalJ2000.Earth_B1(Rho) +
  10:                         VSOP87B_EarthPositionSphericalJ2000.Earth_B2(Rho) + VSOP87B_EarthPositionSphericalJ2000.Earth_B3(Rho) +
  11:                         VSOP87B_EarthPositionSphericalJ2000.Earth_B4(Rho) + VSOP87B_EarthPositionSphericalJ2000.Earth_B5(Rho);
  12:         var radius = VSOP87B_EarthPositionSphericalJ2000.Earth_R0(Rho) + VSOP87B_EarthPositionSphericalJ2000.Earth_R1(Rho) +
  13:                      VSOP87B_EarthPositionSphericalJ2000.Earth_R2(Rho) + VSOP87B_EarthPositionSphericalJ2000.Earth_R3(Rho) +
  14:                      VSOP87B_EarthPositionSphericalJ2000.Earth_R4(Rho) + VSOP87B_EarthPositionSphericalJ2000.Earth_R5(Rho);
  15:         ReferenceCoordinates = new SphericalCoordinates(latitude, longitude, radius);
  16:     };
  17:     Because of = () =>
  18:         ComputedCoordinates = Vsop87OrbitEngine.ComputeSphericalCoordinates(TargetDate,
  19:             SolarSystemBody.Earth,
  20:             ReferenceFrame.EquinoxJ2000);
  21:     It should_match_the_reference_latitude =
  22:         () => ComputedCoordinates.Latitude.ShouldBeCloseTo(ReferenceCoordinates.Latitude, Tolerance);
  23:     It should_match_the_reference_longitude =
  24:         () => ComputedCoordinates.Longitude.ShouldBeCloseTo(ReferenceCoordinates.Longitude, Tolerance);
  25:     It should_match_the_reference_radius =
  26:         () => ComputedCoordinates.Radius.ShouldBeCloseTo(ReferenceCoordinates.Radius, Tolerance);
  27:     static SphericalCoordinates ComputedCoordinates;
  28:     static SphericalCoordinates ReferenceCoordinates;
  29: }
  30:  
  31: [Subject(typeof(Vsop87OrbitEngine), "compute coordinates")]
  32: public class when_computing_rectangular_j2000_coordinates_for_earth : with_target_date_2014_jan_29_midday
  33:     {
  34:     Establish context = () =>
  35:         {
  36:         var x = VSOP87A_EarthPositionRectangularJ2000.Earth_X0(Rho) +
  37:                 VSOP87A_EarthPositionRectangularJ2000.Earth_X1(Rho) +
  38:                 VSOP87A_EarthPositionRectangularJ2000.Earth_X2(Rho) +
  39:                 VSOP87A_EarthPositionRectangularJ2000.Earth_X3(Rho) +
  40:                 VSOP87A_EarthPositionRectangularJ2000.Earth_X4(Rho) +
  41:                 VSOP87A_EarthPositionRectangularJ2000.Earth_X5(Rho);
  42:         var y = VSOP87A_EarthPositionRectangularJ2000.Earth_Y0(Rho) +
  43:                 VSOP87A_EarthPositionRectangularJ2000.Earth_Y1(Rho) +
  44:                 VSOP87A_EarthPositionRectangularJ2000.Earth_Y2(Rho) +
  45:                 VSOP87A_EarthPositionRectangularJ2000.Earth_Y3(Rho) +
  46:                 VSOP87A_EarthPositionRectangularJ2000.Earth_Y4(Rho) +
  47:                 VSOP87A_EarthPositionRectangularJ2000.Earth_Y5(Rho);
  48:         var z = VSOP87A_EarthPositionRectangularJ2000.Earth_Z0(Rho) +
  49:                 VSOP87A_EarthPositionRectangularJ2000.Earth_Z1(Rho) +
  50:                 VSOP87A_EarthPositionRectangularJ2000.Earth_Z2(Rho) +
  51:                 VSOP87A_EarthPositionRectangularJ2000.Earth_Z3(Rho) +
  52:                 VSOP87A_EarthPositionRectangularJ2000.Earth_Z4(Rho) +
  53:                 VSOP87A_EarthPositionRectangularJ2000.Earth_Z5(Rho);
  54:         ReferenceCoordinates = new RectangularCoordinates(x, y, z);
  55:         };
  56:     Because of = () =>
  57:         ComputedCoordinates = Vsop87OrbitEngine.ComputeRectangularCoordinates(TargetDate,
  58:             SolarSystemBody.Earth,
  59:             ReferenceFrame.EquinoxJ2000);
  60:     It should_match_the_reference_x =
  61:         () => ComputedCoordinates.X.ShouldBeCloseTo(ReferenceCoordinates.X, Tolerance);
  62:     It should_match_the_reference_y =
  63:         () => ComputedCoordinates.Y.ShouldBeCloseTo(ReferenceCoordinates.Y, Tolerance);
  64:     It should_match_the_reference_z =
  65:         () => ComputedCoordinates.Z.ShouldBeCloseTo(ReferenceCoordinates.Z, Tolerance);
  66:     static RectangularCoordinates ComputedCoordinates;
  67:     static RectangularCoordinates ReferenceCoordinates;
  68: }
   1: /// <summary>
   2: /// Computes the spherical coordinates for a target body at a given date and with the specified reference frame.
   3: /// </summary>
   4: /// <param name="targetDate">The target date.</param>
   5: /// <param name="solarSystemBody">The solar system body.</param>
   6: /// <param name="referenceFrame">The reference frame.</param>
   7: /// <returns>SphericalCoordinates.</returns>
   8: public static SphericalCoordinates ComputeSphericalCoordinates(double targetDate,
   9:     SolarSystemBody solarSystemBody,
  10:     ReferenceFrame referenceFrame)
  11:     {
  12:     var file = Vsop87DataReader.SelectDataFile(solarSystemBody,
  13:         CoordinateSystem.HeliocentricSphericalCoordinates,
  14:         referenceFrame);
  15:     var vsop87SolutionData = Vsop87DataReader.LoadVsop87DataFromFile(file);
  16:     var latitude = ComputeVsop87Series(targetDate, vsop87SolutionData.VariableData['L']);
  17:     var longitude = ComputeVsop87Series(targetDate, vsop87SolutionData.VariableData['B']);
  18:     var radius = ComputeVsop87Series(targetDate, vsop87SolutionData.VariableData['R']);
  19:     return new SphericalCoordinates(latitude, longitude, radius);
  20:     }
  21:  
  22: /// <summary>
  23: /// Computes the rectangular coordinates for a target body at a given date and with the specified reference frame.
  24: /// </summary>
  25: /// <param name="targetDate">The target date.</param>
  26: /// <param name="solarSystemBody">The solar system body.</param>
  27: /// <param name="referenceFrame">The reference frame.</param>
  28: /// <returns>RectangularCoordinates expressed in Astronomical Units (AU).</returns>
  29:  
  30: public static RectangularCoordinates ComputeRectangularCoordinates(double targetDate,
  31:     SolarSystemBody solarSystemBody,
  32:     ReferenceFrame referenceFrame)
  33:     {
  34:     var file = Vsop87DataReader.SelectDataFile(solarSystemBody,
  35:         CoordinateSystem.HeliocentricRectangularCoordinates,
  36:         referenceFrame);
  37:     var vsop87SolutionData = Vsop87DataReader.LoadVsop87DataFromFile(file);
  38:     var x = ComputeVsop87Series(targetDate, vsop87SolutionData.VariableData['X']);
  39:     var y = ComputeVsop87Series(targetDate, vsop87SolutionData.VariableData['Y']);
  40:     var z = ComputeVsop87Series(targetDate, vsop87SolutionData.VariableData['Z']);
  41:     return new RectangularCoordinates(x, y, z);
  42:     }

…and the coordinate classes…

   1: public struct RectangularCoordinates
   2:     {
   3:     public double X { get; private set; }
   4:     public double Y { get; private set; }
   5:     public double Z { get; private set; }
   6:  
   7:     public RectangularCoordinates(double x, double y, double z) : this()
   8:         {
   9:         X = x;
  10:         Y = y;
  11:         Z = z;
  12:         }
  13:     }
  14:  
  15: public struct SphericalCoordinates
  16:     {
  17:     public SphericalCoordinates(double latitude, double longitude, double radius) : this()
  18:         {
  19:         Latitude = latitude;
  20:         Longitude = longitude;
  21:         Radius = radius;
  22:         }
  23:  
  24:     public double Latitude { get; private set; }
  25:     public double Longitude { get; private set; }
  26:     public double Radius { get; private set; }
  27:     }

Clean Code

The hard coded data that I began with is no longer needed, as it is all read straight from the original data files. In the interests of clean code, I removed it from the solution. Some of the tests aimed to prove that our new implementation, based on this hard coded data, matched the reference implementation. I considered updating those tests to use the data loaded from the data files, but actually what would be the point? We have newer tests that test the end-to-end implementation and verify that it matches the reference implementation. So these intermediate stages are now redundant and I decided to just delete them.

I also took a pass over the entire codebase and sorted out any naming and namespace issues, used ReSharper to ensure consistent formatting and naming and so on. Rather than produce all that code again here, I refer you to my Git server where you can download everything and browse it at your leisure. At this point in time, the latest commit is:
http://stash.teamserver.tigranetworks.co.uk/projects/TNO/repos/ta.orbits/commits/8d3c0c4e7fe717c44509c4ad2514a7f86f009cb5

This is how the final test suite looks, a total of 41 specifications in 11 test contexts:

image

Epilogue

Where next? The solution contains a Markdown file (ReadMe.md) that I have been keeping a ToDo list in. The remaining items in that list are:

  • How to represent the Sun and Earth?
  • Download VSOP87 data from the web if needed (maybe?)
  • Caching strategy, so that each VSOP87 file is loaded only once

In addition, I want the entire OrbitEngine to be ‘pluggable’ so that different algorithms can be used easily, which implies that we need another orbit engine implementation and then to decide on a common interface. Most orbit engines deal with solar system bodies orbiting the Sun (or solar system barycentre) but conceptually I see no reason why there shouldn’t be orbit engines for other solar systems, or for objects orbiting bodies other than the Sun – moons are an obvious case, although I don’t know yet whether an orbit engine would be the right tool for that.

I will look at some of these issues in the next instalment in the series.

Get the Code

All of the code discussed in this article is freely available using anonymous access from our Git server at http://stash.teamserver.tigranetworks.co.uk/projects/TNO/repos/ta.orbits/browse (git clone http://stash.teamserver.tigranetworks.co.uk/scm/tno/ta.orbits.git). Unless otherwise stated, this project is licensed under a Creative Commons Free Culture license

Creative Commons License
Tigra Astronomy Object Oriented Astometry by Tim Long and Tigra Astronomy is licensed under a Creative Commons Attribution 4.0 International License.