Large deviation of latitude values for calculating the point of intersection of arcs on a sphere

Question

My task is to find the intersection point of two arcs on the sphere (if it exists). I used algorithm from http://www.boeing-727.com/Data/fly%20odds/distance.html but in some cases the latitude deviation is too large. What could be causing this?

I have 4 points (lat, lon). I convert the coordinates to Cartesian using the following algorithm:

x = sin(toRadians(p.lat)) * cos(toRadians(p.lon));
y = sin(toRadians(p.lat)) * sin(toRadians(p.lon));
z = cos(toRadians(p.lat));

After that i create the vector V1 (vector director of the first plane, perpendicular to it) and vector V2 (to the second plane)

//coords - class for Cartesian coordinates (x,y,z)
coords V1 = P1 * P2; //Vector multipication (y * rhs.z - rhs.y * z, rhs.x * z - x * rhs.z, x * rhs.y - rhs.x * y)
coords V2 = P3 * P4;

Then i calculate vector director D = V1 * V2; Because the sphere representing the earth has a radius of one, again we calculate the unit vector of D, so that it touches the surface of the sphere. This vector S1 and its opposite S2 directly give the coordinates of the crossing point of the two non-overlapping great circles on the sphere.

length = D.length();
coords S1(D.x / length, D.y / length, D.z / length);
coords S2(-D.x / length, -D.y / length, -D.z / length);

And after that I convert it to spherical coordinates (in Degrees):

lat = toDegrees(atan2(sqrt(x * x + y * y), z));
lon = toDegrees(atan2(y, x));

For example, when crossing the following points (60,30)-(60,60) & (40,50)-(60,50) /*(lat,lon)*/ we get coordinates:

s1: {lat=120.77136585404358 lon=-130.00000000000000 }
s2: {lat=59.228634145956427 lon=50.000000000000014 }

the latitude at the second point is quite different from the correct one (85771.97 meters)

Answer 1

In terms of LATITUDE and LONGITUDE (note that "spherical coordinates" typically use CO-LATITUDE) I think your coordinates should be

x = cos(latitude) * cos(longitude)
y = cos(latitude) * sin(longitude)
z = sin(latitude)

Work these out to get cartesian points P1, P2, P3, P4 on the unit sphere.

Planes through the origin are then defined by their unit normal vectors

V1 = cross(P1,P2)
V2 = cross(P3,P4)

A line of intersection must be perpendicular to both of these, so parallel to cross(V1,V2). As long as this is non-zero you can then normalise this to get a point on the unit sphere.

Finally, convert back to latitude and longitude.

#include <iostream>
#include <cmath>
using namespace std;

const double PI = 4.0 * atan( 1.0 );

//============== Geometry Stuff ========================================

struct Point{ double x, y, z; };

Point operator + ( Point p , Point q  ) { return { p.x + q.x, p.y + q.y, p.z + q.z }; }
Point operator - ( Point p , Point q  ) { return { p.x - q.x, p.y - q.y, p.z - q.z }; }
Point operator * ( double r, Point p  ) { return {   r * p.x,   r * p.y,   r * p.z }; }
Point operator / ( Point p , double r ) { return {   p.x / r,   p.y / r,   p.z / r }; }
double       dot ( Point p , Point q  ) { return p.x * q.x + p.y * q.y + p.z * q.z  ; }
Point      cross ( Point p , Point q  ) { return { p.y * q.z - p.z * q.y, p.z * q.x - p.x * q.z, p.x * q.y - p.y * q.x }; }
double    normsq ( Point p            ) { return dot( p, p ); }
double       abs ( Point p            ) { return sqrt( normsq( p ) ); }

ostream &operator << ( ostream &out, const Point &p ) { return out << p.x << '\t' << p.y << '\t' << p.z; }
istream &operator >> ( istream &in ,       Point &p ) { return in  >> p.x         >> p.y         >> p.z; }

//======================================================================

Point fromLatitudeLongitude( double latitude, double longitude )     // point on unit sphere from latitude, longitude in degrees
{
   double latRad = latitude  * PI / 180.0;
   double lonRad = longitude * PI / 180.0;
   double Clat = cos( latRad ), Slat = sin( latRad );
   double Clon = cos( lonRad ), Slon = sin( lonRad );
   return Point{ Clat * Clon, Clat * Slon, 1.0 * Slat };
}

//======================================================================

void toLatitudeLongitude( Point P, double &latitude, double &longitude )   // latitude, longitude in degrees from point
{
   P = P / abs( P );
   latitude  = asin ( P.z      ) * 180.0 / PI;
   longitude = atan2( P.y, P.x ) * 180.0 / PI;
}

//======================================================================

int main()
{
   double latitude, longitude;
   Point P1 = fromLatitudeLongitude( 60.0, 30.0 );
   Point P2 = fromLatitudeLongitude( 60.0, 60.0 );
   Point P3 = fromLatitudeLongitude( 40.0, 50.0 );
   Point P4 = fromLatitudeLongitude( 60.0, 50.0 );

   Point V1 = cross( P1, P2 );
   Point V2 = cross( P3, P4 );
   Point line = cross( V1, V2 );

   if ( abs( line ) > 1.0e-10 )
   {
      toLatitudeLongitude( line, latitude, longitude );
      cout << "Solution 1: latitude = " << latitude << "      longitude = " << longitude << '\n';

      toLatitudeLongitude( -1.0 * line, latitude, longitude );
      cout << "Solution 2: latitude = " << latitude << "      longitude = " << longitude << '\n';
   }
   else
   {
      cout << "Sorry; unable to solve.\n";
   }
}

For the given example (note: two solutions, on opposite side of the sphere)

Solution 1: latitude = 60.7596      longitude = 50
Solution 2: latitude = -60.7596      longitude = -130

Answer 2

The calculation of spatial coordinates is performed by the formulas:

X=(N+H)*cos(B)*cos(L)
Y=(N+H)*cos(B)*sin(L)
Z=((1-e*e)*N+H)*sin(B)

where

X, Y, Z - spatial rectangular coordinates of the point,

B, L, H -geodetic coordinates of the point,

N - radius of curvature of the first vertical

e - eccentricity of the ellipsoid.

The radius of curvature of the first vertical and the square of the eccentricity ellipsoid are calculated by the formulas:

N=a/sqrt(1-pow(e*sin(B),2))
e = sqrt(2*alpha-pow(alpha,2)

where

a - semi-major axis of the ellipsoid = 1/298,257 84

alpha - ellipsoid contraction = 6 378 136

Data taken from the handbook "ПЗ-90.11"

Answer 3

Your formula for geographical coordinates to xyz does not seem correct. You can use the following code which is meant for ellipsoids, but you can set semi-major axis (a), semi-minor axis (b) to 1.0 you will get your perfect ellipsoid that is a sphere, and if you set height to 0.0 you will be on that sphere's surface:

void geographic_to_xyz(double latitude, double longitude, double height, double a, double b, double &x, double &y, double &z)
{
    double e2_ = (a*a - b*b) / (b*b);
    double c = a*a/b;
    double v1 = sqrt(1 +(e2_ * (cos(latitude) * cos(latitude))));
    double n1 = c / v1;

    x = (n1 + height) * cos(latitude) * cos(longitude);
    y = (n1 + height) * cos(latitude) * sin(longitude);
    z = (pow(b/a,2) * n1 + height) * sin(latitude);
}

Another issue is,even if you use another method, you will most probably have coordinates from WGS84 ellipsoid when using this IRL. Performing spheroidal calculations on ellipsoidal coordinates will yield wrong results. That's why you need to reduce your ellipsoidal latitude first (you only need to reduce latitude because longitude is not affected by the flattening of current ellipsoid), using the formula tan(beta)=(1-f)*tan(phi) where phi is your ellipsoidal latitude, f is flattening of WGS84 ellipsoid (or whichever you use) and beta is the final reduced latitude.