在纬度和经度方面（注意“球坐标”通常使用CO-LATITUDE）我认为你的坐标应该是

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

求出单位球面上的笛卡尔点P1，P2，P3，P4。
然后通过原点的平面由它们的单位法向量定义

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

一条交线必须垂直于这两条线，所以平行于交叉（V1，V2）。只要它不为零，你就可以将它归一化，得到单位球面上的一个点。
最后，转换回纬度和经度。

#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";
   }
}

对于给出的示例（注：两个解决方案，在球体的相对侧）

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

你的xyz地理坐标公式似乎不正确。你可以使用下面的代码，这是椭圆体的意思，但你可以设置半长轴（a），半短轴（b）为1.0，你会得到一个完美的椭圆体，这是一个球体，如果你设置高度为0.0，你将在该球体的表面：

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);
}

另一个问题是，即使您使用其他方法，在使用此IRL时，您也很可能具有来自WGS84椭球的坐标。在椭球坐标上执行椭球计算将产生错误的结果。这就是为什么你需要先减少你的椭球纬度（你只需要减少纬度，因为经度不受当前椭球的展平影响），使用公式tan(beta)=(1-f)*tan(phi)，其中phi是你的椭球纬度，f是WGS84椭球的展平（或任何你使用的），beta是最终减少的纬度。