Finding the first point of great circle Intersection

Finding the first point of great circle Intersection

I have a problem I've been trying to solve and I cannot come up with the
answer. I have written a function in Matlab 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;