Finding the first point of great circle Intersection

Hello Everyone! I have a problem I've been trying to solve and I cannot come up with the answer. I have written a function which, given two lat/lon points and two bearings, will return the two great circle points of intersection.

However, what I really need is the first great circle point of intersection along the two initial headings. I.e. if two airplanes begin at lat/lon points 1 and 2, with initial bearings of bearing1 and bearing2, which of the two great circle intersection points is the first one they encounter? There are many solutions (using haversine) which will give me the closer of the two points, but I don't actually care about which is closer, I care about which I will encounter first given specific start points and headings. There are many cases where the closer of the two intersections is actually the second intersection encountered.

Now, I realize I could do this with lots of conditional statements for handling the different cases, but I figure there's got to be a way to handle it with regard to the order I take all my cross products (function code given below), but I simply can't come up with the right solution! I should also mention that this function is going to be used in a large computationally intensive model, and so I'm thinking the solution to this problem needs to be rather elegant/speedy. Can anyone help me with this problem?

The following is not my function code (I can't list that here), but it is the pseudo code that my function was based off of:

     %Given inputs of lat1,lon1,Bearing1,lat2,lon2,Bearing2:

    %Calculate arbitrary secondary point along same initial bearing from first
    %point
    dAngle = 45;
    lat3 = asind( sind(lat1)*cosd(dAngle) + cosd(lat1)*sind(dAngle)*cosd(Bearing1));
    lon3 = lon1 + atan2( sind(Bearing1)*sind(dAngle)*cosd(lat1), cosd(dAngle)-sind(lat1)*sind(lat3) )*180/pi;

    lat4 = asind( sind(lat2)*cosd(dAngle) + cosd(lat2)*sind(dAngle)*cosd(Bearing2));
    lon4 = lon2 + atan2( sind(Bearing2)*sind(dAngle)*cosd(lat2), cosd(dAngle)-sind(lat2)*sind(lat4) )*180/pi;

    %%Calculate unit vectors
    % We now have two points defining each of the two great circles.  We need
    % to calculate unit vectors from the center of the Earth to each of these
    % points
    [Uvec1(1),Uvec1(2),Uvec1(3)] = sph2cart(lon1*pi/180,lat1*pi/180,1);
    [Uvec2(1),Uvec2(2),Uvec2(3)] = sph2cart(lon2*pi/180,lat2*pi/180,1);
    [Uvec3(1),Uvec3(2),Uvec3(3)] = sph2cart(lon3*pi/180,lat3*pi/180,1);
    [Uvec4(1),Uvec4(2),Uvec4(3)] = sph2cart(lon4*pi/180,lat4*pi/180,1);

    %%Cross product
    %Need to calculate the the "plane normals" for each of the two great
    %circles
    N1 = cross(Uvec1,Uvec3);
    N2 = cross(Uvec2,Uvec4);

    %%Plane of intersecting line
    %With two plane normals, the cross prodcut defines their intersecting line
    L = cross(N1,N2);
    L = L./norm(L);
    L2 = -L;
    L2 = L2./norm(L2);
    %%Convert normalized intersection line to geodetic coordinates

    [lonRad,latRad,~]=cart2sph(L(1),L(2),L(3));
    lonDeg = lonRad*180/pi;
    latDeg = latRad*180/pi;

    [lonRad,latRad,~]=cart2sph(L2(1),L2(2),L2(3));
    lonDeg2 = lonRad*180/pi;
    latDeg2 = latRad*180/pi;

[Cedric]

Actually they might each reach a different point. I might have misunderstood the question, but the way you worded it suggests that both trajectories will converge towards the same point, which is not true. If you imagine the first point being a little after one intersection and the second point being a little after the other intersection, you have a situation were each "plane" will travel towards the intersection that was the closest to the other plane at the beginning.