This answer do all calculations, without boost.
Consider a sphere of radius R = 1.
Points A, B are on a great circle. This great circle gcAB goes also through the center point O of the sphere (required for great circles). Points A, B, O define a plane PL1.
Point P also lies in a great circle.
The minimum distance (measured along the arc of a great circle, not along a 3D straight line) from P to the great circle gcAB is the length of the arc PC.
Plane PL2 of great circle gcPC is perpendicular to the plane PL1.
We want point C, which lies in the line OC, which is the intersection of the two mentioned planes.
.
Plane PL1 is defined by its perpendicular vector pp1. This vector is obtained by the cross product of vectors OA and OB.
Because the plane PL2 is perpendicular to plane PL1, it must contain the vector pp1. So the perpendicular vector pp2 to plane PL2 can be obtained by the cross product of OP and pp1.
A vector ppi in the line OC intersection of both planes is obtained by the cross product of pp1 and pp2.
If we normalize vector ppi and multiply its components by the radius R of Earth, we get the coordinates of point C.
Cross product is not commutative. This means that if we interchange points A,B we get the opposite point C' in the sphere. We can test distances PC and PC' and get their minimum.
To calculate the great-circle distance Wikipedia link for two points A, B, it relies on the angle a between the lines OA and OB.
For best accuracy on all angles we use a = atan2(y, x) where, using radius 1, y= sin(a) and x= cos(a). sin(a) and cos(a) can be calculated by cross product (OA, OB) and dot product (OA, OB) respectively.
Putting all together we have this C++ code:
#include <iostream>
#include <cmath>
const double degToRad = std::acos(-1) / 180;
struct vec3
{
double x, y, z;
vec3(double xd, double yd, double zd) : x(xd), y(yd), z(zd) {}
double length()
{
return std::sqrt(x*x + y*y + z*z);
}
void normalize()
{
double len = length();
x = x / len;
y = y / len;
z = z / len;
}
};
vec3 cross(const vec3& v1, const vec3& v2)
{
return vec3( v1.y * v2.z - v2.y * v1.z,
v1.z * v2.x - v2.z * v1.x,
v1.x * v2.y - v2.x * v1.y );
}
double dot(const vec3& v1, const vec3& v2)
{
return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
}
double GCDistance(const vec3& v1, const vec3& v2, double R)
{
//normalize, so we can pass any vectors
vec3 v1n = v1;
v1n.normalize();
vec3 v2n = v2;
v2n.normalize();
vec3 tmp = cross(v1n, v2n);
//minimum distance may be in one direction or the other
double d1 = std::abs(R * std::atan2(tmp.length() , dot(v1n, v2n)));
double d2 = std::abs(R * std::atan2(tmp.length() , -dot(v1n, v2n)));
return std::min(std::abs(d1), std::abs(d2));
}
int main()
{
//Points A, B, and P
double lon1 = 88.41253929999999 * degToRad;
double lat1 = 22.560206299999997 * degToRad;
double lon2 = 88.36928063300775 * degToRad;
double lat2 = 22.620867969497795 * degToRad;
double lon3 = 88.29580956367181 * degToRad;
double lat3 = 22.71558662052875 * degToRad;
//Let's work with a sphere of R = 1
vec3 OA(std::cos(lat1) * std::cos(lon1), std::cos(lat1) * std::sin(lon1), std::sin(lat1));
vec3 OB(std::cos(lat2) * std::cos(lon2), std::cos(lat2) * std::sin(lon2), std::sin(lat2));
vec3 OP(std::cos(lat3) * std::cos(lon3), std::cos(lat3) * std::sin(lon3), std::sin(lat3));
//plane OAB, defined by its perpendicular vector pp1
vec3 pp1 = cross(OA, OB);
//plane OPC
vec3 pp2 = cross(pp1, OP);
//planes intersection, defined by a line whose vector is ppi
vec3 ppi = cross(pp1, pp2);
ppi.normalize(); //unitary vector
//Radious or Earth
double R = 6371000; //mean value. For more precision, data from a reference ellipsoid is required
std::cout << "Distance AP = " << GCDistance(OA, OP, R) << std::endl;
std::cout << "Distance BP = " << GCDistance(OB, OP, R) << std::endl;
std::cout << "Perpendicular distance (on arc) = " << GCDistance(OP, ppi, R) << std::endl;
}
Wich gives the distances
AP = 21024.4 BP = 12952.1 and PC= 499.493 for the given three points.
Running code here
ta point on the extension of(p, q)which is again an arc because earth is approximately a sphere. It is easy if it is a Cartesian coordinate. But for GPS coordinates its not the case. – Neel Basu