[Proj] Intersection test for two geodesics goes wrong for very long distances

[Mikael Rittri]

Hello all,

I have had a strange problem with geodesics. While I think I have
solved it in a practical sense, by saying: "don't go too near the 
antipode", I wonder exactly how near the antipode it is safe to go,
and why. 
    I post here in the hope that someone would enjoy working on
this kind of problem (hi, Charles). 

The problem was to find the intersection of two radials, where 
a radial is a geodetic line that starts at a given point with 
a given azimuth, but doesn't have a definite end point.

While Carmenta Engine doesn't support this computation directly,
there is support for intersecting two geodetic lines defined by
their start and end points. So, my idea was to compute an arbitrary
end point for each radial, by going forward long enough on the radial.
(This is just the forward or direct geodetic problem.) 

The question then was, how far to go forward? Obviously, one cannot
go longer than 180 degrees or about 20 000 km on the Earth, since
there would then be a shorter geodetic route from start to end on
the backside of the Earth. In fact, since we want to compute on
an ellipsoid, one cannot go quite 180 degrees either. I thought the
maximal safe distance would be  pi * a * (1 - f), where a is the
equatorial radius and f is the flattening. On the WGS84 ellipsoid,
this would be pi * 6378137 * (1 - 1/298.257223563) = 19 970 326 meters. 
(Because, if one starts somewhere on the equator, and goes east
 along the equator to an end point longer than this distance, but 
 shorter than 180 degrees, then the shortest route from start to end
 doesn't follow the equator but takes a shortcut nearer a pole, 
 exploiting the flattening of the ellipsoid.)

So, we tried to let "forward distance" = 19 950 000 meters, and it seemed 
to work for a while. But then we found an example where the intersecting
method wrongly claimed that there was no intersection. It turned out
that this wasn't caused by numerical instability, but by something more
fundamental. 

You see, we have some initial tests in the intersecting method, which
should say when there is definitely no intersection. One condition
is like this, in ASCII graphics:

B1
| .
|   .
|     .
|       . A2             A1
|         -----------------
|       .
|     .
|   .
| .
|
B2 

One geodesic goes from A1 to A2, the other goes from B1 to B2. From
A2, take the azimuths toward A1, B1 and B2. Then if the azimuth
toward A1 doesn't go between the two azimuths toward B1 and B2, 
an intersection would be impossible, I thought. Because A1 and A2
would be on the same side of the geodesic B1B2. 
   Of course, one has to define "between-ness" of azimuths modulo
360 degrees. But the idea is that the "inside" of the angle B1 A2 B2 
would be the side that is less than 180 degrees, and if the line 
from A2 toward A1 doesn't go on the inside, there can be no intersection.

I still feel that this _ought_ to be true. However, here is a counterexample:

The first radial is defined by 
A1 = 17.58139°S 40.17167°W
A1_faz = -90° // A1 forward azimuth

The second radial is defined by 

B1 = 22.23556°S 49.91139°W
B1_faz = 30° // B1 forward azimuth 

They do intersect quite near the start points. 

However, with "forward distance" = 19 950 000 m, the forward computation
gives the end points:

A2 = 17.58126738°N 140.62400233°E

B2 = 22.60088604°N 129.58209305°E

Now, 
Azimuth at A2 toward B1 = 117.78924816°
Azimuth at A2 toward B2 = -62.47659403°

And 117.78924816° - (-62.47659403°) = 180.26584219° which is larger than 180°.
This means that the inside of the angle B1 A2 B2 is south of the angle edges.
But

Azimuth at A2 toward A1 = 90.06652107°

which is north of the A2-toward-B1 azimuth. So, the initial test will say
(wrongly) that there is no intersection.  

In this particular example, I have found that the code finds the intersection
if I lower the forward distance to 19 940 632 m. But there must be some worst
case that requires an even shorter distance, and I don't know how to find
the worst case (except by a lot of trial-and-error). 

I have visualized the counterexample in several map projections and waited
for my brain to say "Aha!", but in vain. I think the intuition behind the
initial test is that the line A1A2 is on one (oblique) hemisphere, and the
triangle B1B2A2 is on the other hemisphere. It feels right, but clearly it
isn't true. 

The pragmatic solution seems to be to use a much shorter forward distance,
something like 15 000 000 m perhaps, but I would really like to understand
what's going on here. 

Credit: the values above have been computed on the GeographicLib site,
http://geographiclib.sourceforge.net/cgi-bin/Geod

Regards,

Mikael Rittri
Carmenta
Sweden
http://www.carmenta.com