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

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"