Just another WordPress.com site

WGS84 Ellipsoid Calculations

(source code links below)
12/7/2010, added support for building as a DLL and using within Microsoft Excel VBA.

In late 2008, I was in negotiations with an aviation company to write some navigation design software. I use to work for the company and was already familiar with their software stacks and engineering processes. Quite a few of the applications the company uses, I had written or helped implement. The chances of actually getting the contract were slim though. Internal growth of the company was fast, and different departments were motivated to go a different direction than the department I was in talks with.

Even though the chances of getting the contract were slim, I wanted to get a jump on some crucial code that would be required from the onset. So I set out on October 15, 2008, to write all the formulas described in FAA Order 8260.54A, Appendix 2, and a test harness for the data published in attachment C of the same document.

Initially I figured it would take 2 to 4 weeks to write the library and test harness based on the completeness of 8260.54A. Almost immediately, I discovered the document to have many errors and omissions, which caused days long investigations into why the formulas didn’t produce the expected test results (the projects code, comments and file Application test Results.txt document some of the problems and how they are resolved).

On December 18, 2008, I finished the work, had no contract, but was proud that I’d resolved each and every formula in Appendix 2 and they produced the expected test results. Now, a year and a half later, I’m releasing the code for the formulas and test harness.

I wrote a standalone WGS84 Geodesic Calculator (beta) for Windows based on the formulas, and some people from the GeoCache community have used it to create and solve complicated puzzles.

The hardest part about releasing the code was which license to choose. Ultimately I decided for the liberal Apache License, Version 2.0, and it’s my contribution to the open source community.

WGS84 Ellipsoid Calculation Types

Arc Intercept Geodesic Locus Intersect Point on LocusP
Course Intersect Geometric Mean Meridional Point to Arc Tangents
Direct Arc Length Locus Arc Intersect Pt is on Arc
Discretized Arc Length Locus Course at Point Pt is on Geodesic
Distance Intersection Locus Intersect Pt is on Locus
Distance to LocusP Locus Perpendicular Intercept Tangent Fixed Radius Arc
Distance to LocusD Locus Tangent Fixed Radius Arc Vincenty Direct
ECEF Perpendicular Intercept Vincenty Inverse
Geodesic Arc Intercept Perpendicular Tangent Points

There are a few other implemented formulas which are support to the ones listed above.  Also, the TerpsTest code for validating the formulas is a good source to see how to use the formulas in your own code.

Available Source Code

The GeoFormulas and TerpsTest (test harness) code is hosted on GitHub, and released under the Apache License, Version 2.0. Note: TerpsTest requires Boost::RegEx. Outside of STL, GeoFormulas is not dependent on external libraries.

The source code for the standalone Windows WGS84 Geodesic Calculator is not available at this time.


July 14, 2010 Posted by | Ellipsoid, FAA, WGS84 | , , , , , , , , | 3 Comments

N-dimension array initialization

In the previous article Container assignment list I show how to assign items to STL containers using a list of items separated by a comma.

This example will show how to do the same with a N-dimensional array, with debug build bounds checking. First an example which shows how to declare and initialize a N-dimensional array of type int or type double with a comma separated list of items.

  1. // define a 3 dimensional array of int's
  2. // and initialize.
  3. MyArray<int, 3> intArray3D;
  4. intArray3D = 1, 2, 3;
  6. // define a 2 diemsional array of double's
  7. // and initialize.
  8. MyArray<double, 2> doubleArray2D;
  9. doubleArray2D = 10.0, 20.0;
  11. // define a 3 x 3 matrix of double and initialize
  12. MyArray<double, 9> doubleMatrix3x3;
  13. doubleMatrix3x3 = 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0;

The N-Dimensional array class implementation

The array class code below only shows enough code to get the assignment and comma operator overloaded. Other useful code and overloads that this class would have is not shown, in order to keep the discussion simple and on topic.

MyArray Class
  1. template<class T, int ArraySize>
  2. class MyArray
  3. {
  4. public:
  5.     MyArray()
  6.     {
  7.         assert(ArraySize != 0);
  8.         // Set items in array to 0
  9.         for(size_t i = 0; i < ArraySize; ++i)
  10.             m_array[i] = 0;
  11.     }
  13.     Initializer<T, ArraySize> operator=(T val)
  14.     {
  15.         // store RHS of = to first array item location
  16.         m_array[0] = val;
  18.         // Call Initializer construct with address of m_array[1],
  19.         // return Initializer object.
  20.         // Allows the overloaded comma operator for Initializer
  21.         // to be called repeatedly.
  22.         return Initializer<T, ArraySize>(&m_array[1]);
  23.     }
  24. private:
  25.     T m_array[ArraySize];
  26. };

MyArray is a template class that needs 2 template parameters, the first is the type of array to create, and the second is the number of elements the array can store.

The default constructor assigns 0 to each element in the array, if this does not make sense for the type object being stored in the array then modify as needed.

The overloaded assignment operator takes the first item in the list, right of the = sign (RHS), and assigns it to first position in the array. It then calls Initializer’s constructor with the address of the second element in the array, and returns the Initializer object.

This sets up Initializer, which overloads the comma operator, to handle the rest of the assignment list. Overloading the assignment and comma operators is discussed in the article Container assignment list.

The Initializer class is simple and documented below.

Initializer Class for MyArray
  1. // Initializer class. Helper class for MyArray,
  2. // keeping track of the current element being
  3. // assigned, with the overloaded comma operator,
  4. template<class T, int ArraySize>
  5. class Initializer
  6. {
  7. public:
  8.     Initializer() : m_ptr(0)
  9.     {
  10.         #ifndef NDEBUG
  11.             m_nCount = 0;
  12.         #endif
  13.     }
  15.     // Constructor called by MyArray operator=
  16.     // ptr will equal MyArray::m_array[1]. Since
  17.     // For debug builds set m_nCount to 1 since
  18.     // MyArray already assigned element [0].
  19.     Initializer(T * ptr) : m_ptr(ptr)
  20.     {
  21.         #ifndef NDEBUG
  22.             m_nCount = 1;
  23.         #endif
  24.     }
  26.     // Overloaded comma operator
  27.     Initializer & operator,(T val)
  28.     {
  29.         #ifndef NDEBUG
  30.             // Debug assert for array overrun
  31.             assert(m_nCount < ArraySize);
  32.             ++m_nCount;
  33.         #endif
  35.         // store value in pointer to array, then
  36.         // increment pointer to array to next location
  37.         *m_ptr++ = val;
  38.         return *this;
  39.     }
  40. private:
  41.     // Pointer to array location
  42.     T * m_ptr;
  43.     #ifndef NDEBUG
  44.         int m_nCount;
  45.     #endif
  46. };

July 12, 2010 Posted by | Arrays, C++ Programming, Templates | , , , , , | Leave a comment

Container assignment lists

A few weeks ago while reading the Blitz++ Documentation I noticed Blitz++ was initializing Matrices and other containers using commas.

  1. Container<int> bucket;
  2. bucket = 1, 2, 3, 4, 5;

Finding this both useful and cool, had to figure out how it was done.

But first…

If the container has a std::vector member then a constructor taking both the start and end of the array can be used to initialize vector. First an array is initialized, and the beginning and end of the array is passed to the container’s constructor.

  1. const int elements = 5;
  2. int a[elements] = {1, 2, 3, 4, 5, };
  3. Container<int> vec(a, &a[elements]);

Although this is doable it isn’t flexible and you might be incline to code a bunch of push_back calls

  1. Container<int> vec;
  2. vec.push_back(1);
  3. vec.push_back(2);
  4. vec.push_back(3);
  5. vec.push_back(4);
  6. vec.push_back(5);

If we have to go that route it would be much nicer to write code as shown at the beginning of the article.

Overloading a couple operators

The goal again, is to be able to write code such as

  1. Container<int> bucket;
  2. bucket = 1, 2, 3, 4, 5;

which assigns values to a container using a list.

Below is code for a template container class that has std::vector as a member. Only the two operators to be overloaded are shown in the example.

Container class
  1. template<class T>
  2. class Container
  3. {
  4. public:
  5. std::vector<T> & operator=(const T & nVal) {
  6. m_vector.clear();
  7. m_vector.push_back(nVal);
  8. return m_vector;
  9. }
  10. friend std::vector<T> & operator,(std::vector<T> & lst, T val)
  11. {
  12. lst.push_back(val);
  13. return lst;
  14. }
  15. std::vector<T> m_vector;
  16. };

Assignment operator

The assignment operator is overloaded to take a const reference T as the first parameter and returns a reference to m_vector. In the previous examples T is an int.  When the compiler sees

  1. bucket = 1

it grabs the value to the right of = and if the correct type, produces code for the assignment operator. The overloaded operator= first clears the contents of m_vector, then pushes the value onto m_vector.

The assignment operator returns a reference of m_vector which allows the overloaded comma operator to be called from the assignment list.

Comma operator

The second operator overloaded for the container class is the comma operator, which evaluates the object type to the left, evaluates the object type to the right and operates on it.

For the purposes of the example, the assignment operator returned a reference to std::vector<int> (m_vector) from the first item in the list.

  1. bucket = 1, 2, 3, 4, 5;

The compiler uses the returned type on the first seen comma (LHS), looks to the right of the comma and sees it is an int type (RHS). This matches the overloaded comma operator’s definition of

  1. friend std::vector<int> & operator,(std::vector<int> & lst, int val)

Since the comma operator is a binary operator, it needs have global scope and declaring the function as a friend provides that.

For the container class the overloaded comma operator simply pushes the RHS hand value onto the LHS object, and returns a reference to the std::vector<T> (m_vector), and repeats until no items in the list are left.


This type of operator overloading can eliminate ugly lists of repeated push_back calls and can easily be used with the other STL containers.

It is important that your container class has a member of the STL container you want to use. Do not inherit from one of the STL containers. They do not have a virtual destructor, which is needed for a derived class.

July 9, 2010 Posted by | C++ Programming, STL Containers, Templates | , , , , , | 1 Comment