Feature/geographic

The purpose of this PR is to move an algorithm resolving vincenty inverse problem outside of the vincenty strategy.

The goal was to add the algorithm a way of calculating not only distance but also forward and reverse azimuths and in the same time to reuse as much code as possible. I implemented a class performing a common part of the algorithm in the ctor and allowing to calculate desired result with member functions. This way it's possible to e.g. calculate the distance and azimuth without recalculating most of the intermediate values:

    vincenty_inverse<CT> vi(lon1, lat1, lon2, lat2, spheroid);
    vi.distance();
    vi.azimuth12();

The next step is to move the vincenty and andoyer strategies out of the extensions. Since Vincenty's algorithm could be used one way or another in various strategies would it be a good idea to rename it (e.g. to `vincenty_distance`)? Or would you prefer to leave the name `vincenty` for a distance strategy and pick some suitable names for possible other strategies in the future?
AFAIU it's said that Vincenty's formula allows to calculate the distance with a millimeter accuracy. Sigma is an angle in radians and it's used later in calculations so I'm assuming it should be calculated with **at least** millimeter precision. On the equator (worst case scenario) 1e-12 more or less corresponds to 0.006mm (1mm~1e-10). I'm guessing that this threshold is some kind of empirical choice. We shouldn't choose value too small to not unnecessarily increase the number of iterations. So it's a precision vs performance tradeof.

With that said I think that the true reason why this threshold was chosen is that it's mentioned on one of the pages in the comment (http://www.movable-type.co.uk/scripts/latlong-vincenty.html). On another one (http://futureboy.homeip.net/fsp/colorize.fsp?fileName=navigation.frink) 1e-6 arcsec is used which is more or less equal to 5e-12 rad if I'm correct. Btw, PostGIS implementation of Vincenty uses 1e-9.

Yes, it'd be reasonable to use a threshold suitable for a floating point type used for calculation. I'm not sure if we should decrease it to achieve more precision but for sure we could consider using greater threshold for floats since the machine epsilon in this case is 10k times greater (~1e-8), so indeed we could use epsilon here for optimization purposes. Another solution would be to always use at least double and then cast the result back, as it's done for integral types.