Opened 6 years ago

Closed 2 years ago

# command-line geod fails when listing points along the equator

Reported by: Owned by: mikebanahan warmerdam major CommandlineUtilities unspecified geod, geodesic

### Description

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.

### comment:1 Changed 3 years ago by karney

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

### comment:2 Changed 2 years ago by warmerdam

• Resolution set to fixed
• Status changed from new to closed

The code from #197 (reimplementing geod) is now checked in and the reported request gives the following which I suppose is appropriate:

```src/geod +ellps=WGS84 +n_S=25 +lat_1=0dN +lon_1=-2dW +lat_2=0dN +lon_2=2dE -f %f
0.000000	-2.000000
0.000000	-1.840000
0.000000	-1.680000
0.000000	-1.520000
0.000000	-1.360000
0.000000	-1.200000
0.000000	-1.040000
0.000000	-0.880000
0.000000	-0.720000
0.000000	-0.560000
0.000000	-0.400000
0.000000	-0.240000
0.000000	-0.080000
0.000000	0.080000
0.000000	0.240000
0.000000	0.400000
0.000000	0.560000
0.000000	0.720000
0.000000	0.880000
0.000000	1.040000
0.000000	1.200000
0.000000	1.360000
0.000000	1.520000
0.000000	1.680000
0.000000	1.840000
0.000000	2.000000
```
Note: See TracTickets for help on using tickets.