Этот ответ делают все вычисления, без повышения .
Рассмотрим сферу радиуса R = 1.
Точки A, B находятся на большом круге . Этот большой круг gcAB проходит также через центральную точку O сферы (требуется для больших кругов). Точки A , B , O определяют плоскость PL1.
Точка P также лежит в большом круге.
Минимальное расстояние (измеренное по дуге большого круга, а не по прямой 3D) от P до большого круга gcAB - это длина дуги PC .
Плоскость PL2 большого круга gcPC перпендикулярна плоскости PL1 .
Нам нужна точка C , которая лежит на линии OC , которая является пересечением двух упомянутых плоскостей.
.
Плоскость PL1 определяется ее перпендикулярным вектором pp1. Этот вектор получается из перекрестного произведения векторов OA и OB.
Поскольку плоскость PL2 перпендикулярна плоскости PL1 , она должна содержать вектор pp1. Таким образом, перпендикулярный вектор pp2 к плоскости PL2 может быть получен путем перекрестного произведения OP и pp1.
Вектор ppi в линии OC пересечения обеих плоскостей получается путем перекрестного произведения pp1 и pp2.
Если мы нормализуем вектор ppi и умножим его составляющие на радиус R Земли, мы получим координаты точки C .
Кросс-произведение не является коммутативным. Это означает, что если мы поменяем точки A, B, мы получим противоположную точку C ' в сфере. Мы можем проверить расстояния PC и PC' и получить их минимум.
Для расчета расстояния по большому кругу Ссылка на Википедию для двух точек A , B , она опирается на угол a между строк OA и OB.
Для лучшей точности на всех углах мы используем a = atan2(y, x), где, используя радиус 1, y= sin(a) и x= cos(a). sin(a) и cos(a) можно рассчитать по перекрестному произведению (OA, OB) и точечному произведению (OA, OB) соответственно.
Собрав все вместе, мы имеем этот код C ++:
#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;
}
который дает расстояния
AP = 21024,4 BP = 12952,1 и PC = 499,493 для данных трех точек.
Рабочий код здесь