找到使用球形坐标表达的线段之间的距离

我试图找到以度（纬度，经度）表示的两个linegements之间的距离。例如，

LINESTRING1 is ((4.42553, -124.159), (4.42553, -124.159)) and
LINESTRING2 is ((-0.061225, -128.428), (-0.059107, -128.428))

然而，由于地球的球形尺寸，找到此距离的通常方法使我距离不正确。是否有某种方式可以改进下面的代码，以使我返回正确的距离，即：558 km

#include <cstdlib>
#include <map>
#include <vector>
#include <algorithm>
#include <string>
#include <future>
#include <string>
#include <fstream>
#include <sstream>
#include <iostream>
using namespace std;
double dot( double* V1, double* V2 )
{
    return V1[0]*V2[0] + V1[1]*V2[1];
}
void lineLineDist2()
        //Compute the square of the minimum distance between two line segments
{
        double seg1[4];
        seg1[0]=4.42553; seg1[1]=-124.159; seg1[2]=4.42553; seg1[3]=-124.159;
        double seg2[4];
        seg2[0]=-0.061225; seg2[1]=-128.428; seg2[2]=-0.059107; seg2[3]=-128.428;

        double u[2],v[2],w[2];
        //Vector   u = S1.P1 - S1.P0;
        u[0] = seg1[2] - seg1[0];
        u[1] = seg1[3] - seg1[1];
        //Vector   v = S2.P1 - S2.P0;
        v[0] = seg2[2] - seg2[0];
        v[1] = seg2[3] - seg2[1];
        //Vector   w = S1.P0 - S2.P0;
        w[0] = seg1[0] - seg2[0];
        w[1] = seg1[1] - seg2[1];
        double    a = dot(u,u);         // always >= 0
        double    b = dot(u,v);
        double    c = dot(v,v);         // always >= 0
        double    d = dot(u,w);
        double    e = dot(v,w);
        double    D = a*c - b*b;        // always >= 0
        double    sc, sN, sD = D;       // sc = sN / sD, default sD = D >= 0
        double    tc, tN, tD = D;       // tc = tN / tD, default tD = D >= 0                
        // compute the line parameters of the two closest points
        if (D < 0.00000001) {           // the lines are almost parallel
                sN = 0.0;               // force using point P0 on segment S1
                sD = 1.0;               // to prevent possible division by 0.0 later
                tN = e;
                tD = c;
        }
        else {                          // get the closest points on the infinite lines
                sN = (b*e - c*d);
                tN = (a*e - b*d);
                if (sN < 0.0) {         // sc < 0 => the s=0 edge is visible
                        sN = 0.0;
                        tN = e;
                        tD = c;
                }
                else if (sN > sD) {     // sc > 1  => the s=1 edge is visible
                        sN = sD;
                        tN = e + b;
                        tD = c;
                }
        }
        if (tN < 0.0) {                 // tc < 0 => the t=0 edge is visible
                tN = 0.0;
                // recompute sc for this edge
                if (-d < 0.0)
                        sN = 0.0;
                else if (-d > a)
                        sN = sD;
                else {
                        sN = -d;
                        sD = a;
                }
        }
        else if (tN > tD) {             // tc > 1  => the t=1 edge is visible
                tN = tD;
                // recompute sc for this edge
                if ((-d + b) < 0.0)
                        sN = 0;
                else if ((-d + b) > a)
                    sN = sD;
                else {
                        sN = (-d +  b);
                        sD = a;
                }
        }
        // finally do the division to get sc and tc
        sc = (abs(sN) < 0.00000001 ? 0.0 : sN / sD);
        tc = (abs(tN) < 0.00000001 ? 0.0 : tN / tD);
        // get the difference of the two closest points
        double dp[2];
        //Vector   dP = w + (sc * u) - (tc * v);  // =  S1(sc) - S2(tc)
        dp[0] = w[0] + sc *u[0] - tc*v[0];
        dp[1] = w[1] + sc *u[1] - tc*v[1];                
        cout<<"distance="<<dot(dp,dp)<<"n";
}
//---------------------------------------------------------------------------
int main()
{
    lineLineDist2();
    return 0; 
}