As of Version 7.0, geodesy functionality is built into the Wolfram Language.

Geodesy Package

Geodesy is the branch of science that deals with such topics as determining positions and areas over large parts of the Earth, the shape and size of the Earth, and the variations in the Earths gravitational and magnetic fields. The primary functions in this package are used for determining the distance between two points.

SphericalDistance[pos1,pos2]calculate the distance between the posi using a spherical model
SpheroidalDistance[pos1,pos2]calculate the distance between the posi using a spheroidal model

Finding the distance between two points on a sphere.

Each position can be given in degrees as a latitudelongitude pair. A coordinate can also be expressed in the form {degrees,minutes,seconds}. A negative value for a coordinate indicates that the coordinate is south latitude or west longitude. Distances are returned in kilometers as the default.

This loads the package.
Click for copyable input
This gives the distance between Null lat., Null long., and Null N. lat., Null E. long., using a spherical model.
Click for copyable input
Here is the distance between the same points using a spheroidal model.
Click for copyable input
This is the difference between the two models.
Click for copyable input

The spherical model is based on an exact formula, while the spheroidal model is an approximation algorithm that is fairly good for distances of less than 10,000 kilometers. Because of the nature of this approximation, all computation with SpheroidalDistance is done with machineprecision numbers.

Radius->radiusspecify the radius of the sphere in the spherical model
SemimajorAxis->lengthspecify the length of the semimajor axis in the spheroidal model
Eccentricity->valuespecify the value of the eccentricity in the spheroidal model

Options to control the size of the shape in the models.

The distances between points can be modified by using different values for the radius of the sphere, or the length of the semimajor axis and the eccentricity of the spheroid using the options Radius, SemimajorAxis, and Eccentricity, respectively. The default values for these options are those for Earth from the WGS84 standard, in kilometers.

ToAuthalicRadius[semimajor,eccentricity]compute radius of authalic sphere given spheroids semimajor axis and eccentricity
GeodeticToAuthalic[{lat,long},eccentricity]compute coordinates of corresponding point on authalic sphere given latitude and longitude of point on spheroid with specified eccentricity
ToDegrees[{deg,min,sec}]convert coordinate in degreeminutesecond format to degrees
ToDMS[deg]convert coordinate in degrees to degreeminutesecond format

Conversions between models and coordinate systems.

The simplest way to convert to the sphere from the spheroid is to use the authalic sphere. This sphere has the same surface area as the reference spheroid.

Here is a comparison of the two shapes, with an ellipse whose eccentricity is much larger than the Earth, for effect.
Click for copyable input

Note that the eccentricity used in the diagram was .6. The actual eccentricity of Earths spheroid is approximately .081.