command-line geod fails when listing points along the equator

[mikebanahan]

Using geod to list points on a great circle route around the equator causes NaN values to be produced, e.g

geod +ellps=WGS84 +n_S=25 +lat_1=0dN +lon_1=-2dW +lat_2=0dN +lon_2=2dE -f %f
produces
0.000000 -2.000000
nan nan
nan nan
... etc

Looking at the code there are two clear problems. The first is not relevant to this but should be fixed. Line 18 in geod_for.c reads

if ((merid = fabs(sina12 = sin(al12)) < MERI_TOL)) {

whereas I believe it SHOULD be

if ((merid = (fabs(sina12 = sin(al12)) < MERI_TOL))) {

i.e instead of assigning the fabs value to merid, it should assign the truth value.

That then prompts what appears to be correct behaviour on north-south geodesics along the 0 longitude meridian

The cause of NaN is at lines 41-46:

41 if (merid) s1 = HALFPI - th1;
42 else {
43 s1 = (fabs(M) >= 1.) ? 0. : acos(M);
44 s1 = sinth1 / sin(s1);
45 s1 = (fabs(s1) >= 1.) ? 0. : acos(s1);
46 }

Where on line 43 s1 becomes zero and line 44 immediately produces divide-by-zero

As a result I've had to switch to using Geographic Library for great-circle work, which is a shame. I don't understand the code well enough (too few comments here and what are to me, opaque variable names) to fix it. It doesn't help that I'm not a good enough geographer/mathematician to understand the algorithm of course.

[karney]

Note that the patch supplied with ticket #197 fixes this problem.