title: Find distance between line segments expressed using spherical coordinates
body: <p>I am trying to find Distance between two linesegments expressed in degrees(latitude, longitude). For example,</p>
  
  <pre><code>LINESTRING1 is ((4.42553, -124.159), (4.42553, -124.159)) and
  LINESTRING2 is ((-0.061225, -128.428), (-0.059107, -128.428))
  </code></pre>
  
  <p>However, the usual way to find this distance is giving me incorrect distance because of the spherical dimensions of the earth. Is there some way by which I may improve the code below such that it returns me the correct distance which is: 558 km </p>
  
  <pre><code>#include &lt;cstdlib&gt;
  #include &lt;map&gt;
  #include &lt;vector&gt;
  #include &lt;algorithm&gt;
  #include &lt;string&gt;
  #include &lt;future&gt;
  #include &lt;string&gt;
  #include &lt;fstream&gt;
  #include &lt;sstream&gt;
  #include &lt;iostream&gt;
  
  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 &gt;= 0
          double    b = dot(u,v);
          double    c = dot(v,v);         // always &gt;= 0
          double    d = dot(u,w);
          double    e = dot(v,w);
          double    D = a*c - b*b;        // always &gt;= 0
          double    sc, sN, sD = D;       // sc = sN / sD, default sD = D &gt;= 0
          double    tc, tN, tD = D;       // tc = tN / tD, default tD = D &gt;= 0                
  
          // compute the line parameters of the two closest points
          if (D &lt; 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 &lt; 0.0) {         // sc &lt; 0 =&gt; the s=0 edge is visible
                          sN = 0.0;
                          tN = e;
                          tD = c;
                  }
                  else if (sN &gt; sD) {     // sc &gt; 1  =&gt; the s=1 edge is visible
                          sN = sD;
                          tN = e + b;
                          tD = c;
                  }
          }
  
          if (tN &lt; 0.0) {                 // tc &lt; 0 =&gt; the t=0 edge is visible
  
                  tN = 0.0;
  
                  // recompute sc for this edge
                  if (-d &lt; 0.0)
                          sN = 0.0;
                  else if (-d &gt; a)
                          sN = sD;
                  else {
                          sN = -d;
                          sD = a;
                  }
          }
          else if (tN &gt; tD) {             // tc &gt; 1  =&gt; the t=1 edge is visible
  
                  tN = tD;
  
                  // recompute sc for this edge
                  if ((-d + b) &lt; 0.0)
                          sN = 0;
                  else if ((-d + b) &gt; a)
                      sN = sD;
                  else {
                          sN = (-d +  b);
                          sD = a;
                  }
          }
  
          // finally do the division to get sc and tc
          sc = (abs(sN) &lt; 0.00000001 ? 0.0 : sN / sD);
          tc = (abs(tN) &lt; 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&lt;&lt;"distance="&lt;&lt;dot(dp,dp)&lt;&lt;"\n";
  }
  //---------------------------------------------------------------------------
  int main()
  {
      lineLineDist2();
      return 0; 
  }
  </code></pre>