Using a spherical model of the earth may give adequate accuracy and leads to simple fast calculations.

Convert all coordinates into earth-centered (3D) cartesian coordinates. For example, the formula

(cos(lon)*cos(lat), sin(lon)*cos(lat), sin(lat))

will do. (It uses a distance measure in which the earth's radius is one unit, which is convenient.)

Writing X0 = (x0,y0,z0) for the start point and X1 = (x1,y1,z1) for the destination point, which define the great circle (provided X0 is distinct from X1 and the two are not diametrically opposite), let U be the normalized cross product of X0 and X1. This is computed in two steps:

V = (xv, yv, zv) = (y0*z1 - z0*y1, z0*x1 - x0*z1, x0*y1 - y0*x1)

The length of V is

|V| = sqrt(xv^2 + yv^2 + zv^2)

Normalization stretches V to unit length:

U = (xu, yu, zu) = V / |V| = (xv/|V|, yv/|V|, zv/|V|).

The oriented 3D distance between any point X = (x,y,z) and the plane of this great circle is just the dot product of X with Z, given by

d = X * U = x*xu + y*yu + z*zu

The closest point in terms of the distance on the earth's surface is the one that is closest to the plane: thus, it has the smallest absolute value of d.

This figure shows a great circle (in black) determined by the two white points and 2000 random points on the sphere colored and shaded according to their absolute 3D distance to the plane of that great circle; that is, |d|.

Having found a closest point, project it to the great circle by first projecting it to the great circle's plane (in 3D) and then extending that radially outward to the earth's surface. The projection merely subtracts d * U:

X' = (x', y', z') = X - d*U = (x - d*xu, y - d*yu, z - d*zu).

The radial projection simply renormalizes X' in the same way V was renormalized to U:

X'' = X' / |X'|.

(This will be problematic if |X'| = 0, which happens when the closest point is one of the poles of the great circle. Include a test in the code for this condition, if it could happen, and deal with it separately, using the sign of d to identify which pole.)

If desired, convert the coordinates of X'' back to (lat, lon) using the usual formulas.

Using a spherical model of the earth may give adequate accuracy and leads to simple fast calculations.

Convert all coordinates into earth-centered (3D) cartesian coordinates. For example, the formula

(cos(lon)*cos(lat), sin(lon)*cos(lat), sin(lat))

will do. (It uses a distance measure in which the earth's radius is one unit, which is convenient.)

Writing X0 = (x0,y0,z0) for the start point and X1 = (x1,y1,z1) for the destination point, which define the great circle (provided X0 is distinct from X1 and the two are not diametrically opposite), let U be the normalized cross product of X0 and X1. This is computed in two steps:

V = (xv, yv, zv) = (y0*z1 - z0*y1, z0*x1 - x0*z1, x0*y1 - y0*x1)

The length of V is

|V| = sqrt(xv^2 + yv^2 + zv^2)

Normalization stretches V to unit length:

U = (xu, yu, zu) = V / |V| = (xv/|V|, yv/|V|, zv/|V|).

The oriented 3D distance between any point X = (x,y,z) and the plane of this great circle is just the dot product of X with Z, given by

d = X * U = x*xu + y*yu + z*zu

The closest point in terms of the distance on the earth's surface is the one that is closest to the plane: thus, it has the smallest absolute value of d.

This figure shows a great circle (in black) determined by the two white points and 2000 random points on the sphere colored and shaded according to their absolute 3D distance to the plane of that great circle; that is, |d|.

Having found a closest point, project it to the great circle by first projecting it to the great circle's plane (in 3D) and then extending that radially outward to the earth's surface. The projection merely subtracts d * U:

X' = (x', y', z') = X - d*U = (x - d*xu, y - d*yu, z - d*zu).

The radial projection simply renormalizes X' in the same way V was renormalized to U:

X'' = X' / |X'|.

(This will be problematic if |X'| = 0, which happens when the closest point is one of the poles of the great circle. Include a test in the code for this condition, if it could happen, and deal with it separately, using the sign of d to identify which pole.)

If desired, convert the coordinates of X'' back to (lat, lon) using the usual formulas.