Calculating intersection of geodesics on ellipsoid (WGS84, ПЗ-90)

Question

I found the formulas for finding the circle point on the sphere on the site http://www.boeing-727.com/Data/fly%20odds/distance.html, but I can’t figure out how to transfer the calculations to the ellipsoid (for example, wgs84 or ПЗ-90). How big is the error at latitudes up to 80 degrees, if you use a spherical model of the earth. And if the discrepancies are large (more than a meter), then how to find the intersection points, taking into account the parameters of the earth?

I tried to use the formulas for obtaining Cartesian coordinates, given in the PZ-90 reference book, but together with the formula on the site given, they give the wrong result.

Algorithm from http://www.boeing-727.com/Data/fly%20odds/distance.html

LatLonSpherical LatLon_intersection(LatLonSpherical p1, LatLonSpherical p2, LatLonSpherical p3, LatLonSpherical p4) {
    //math from here: http://www.boeing-727.com/Data/fly%20odds/distance.html#arc2
    coords P1(p1), P2(p2), P3(p3), P4(p4);//Преобразование в декартову систему координат
    coords V1 = cross(P1, P2); //Vector multipication
    coords V2 = cross(P3, P4);
    //Calculate unit vector for V1 V2
    double length = V1.length();
    coords U1(V1.x / length, V1.y / length, V1.z / length);
    length = V2.length();
    coords U2(V2.x / length, V2.y / length, V2.z / length);
    //vector director
    //coords D = U1 * U2;
    coords D = cross(V1, V2);
    length = D.length();
    //Vector S1 and its opposite S2 directly give the coordinates of the crossing point of the two non-overlapping great circles on the sphere.
    coords S1(D.x / length, D.y / length, D.z / length);
    coords S2(-D.x / length, -D.y / length, -D.z / length);
    //Converting Cartesian system into geographical location    
    LatLonSpherical point1(D.x, D.y, D.z);
    //Check location of crossing points with arcs of great circles
    LatLonSpherical s1(S1.x, S1.y, S1.z), s2(S2.x, S2.y, S2.z);
    if (abs(p1.distanceTo(p2) - p1.distanceTo(s1) - p2.distanceTo(s1)) < 1 &&   
        abs(p3.distanceTo(p4) - p3.distanceTo(s1) - p4.distanceTo(s1)) < 1) {
        return s1; //s1 on arc
    }
    if (abs(p1.distanceTo(p2) - p1.distanceTo(s2) - p2.distanceTo(s2)) < 1 &&
        abs(p3.distanceTo(p4) - p3.distanceTo(s2) - p4.distanceTo(s2)) < 1) {
        return s2; //s2 on arc
    }
    s2.error = -1;//Возврат ошибки
    return s2;
}

Answer 1

FAA document (https://www.faa.gov/documentLibrary/media/Order/8260.54A.pdf) on page A2-32 is algorithm for computing intersection point of two geodesics