Version 1.46 returns NaN, but old version 1.45 works fine

[zed]

It's all about C version of library.

I have simple test code, that works fine with v.1.45, but it didn't work with current release:

int GeodSolve74() {
  double lat2, lon2, azi2;
  struct geod_geodesic g;
  int result = 0;
  geod_init(&g, wgs84_a, 0);
  //geod_init(&g, wgs84_a, wgs84_f);
  geod_direct(&g, -3.688856, -104.23828, 0, 4000000,
              &lat2, &lon2, &azi2);  
  result += assertEquals(lat2, 32.243755, 0.5e-5);
  result += assertEquals(lon2, -104.23828, 0.5e-5);
  result += assertEquals(azi2, 0.0, 0.5e-5);
  return result;
}

It works on WGS84 ellipsoid, but I want to calc coordinates on Google Sphere (Radius A = Radius B = 6378137 http://spatialreference.org/ref/sr-org/7483/).

[Charles Karney]

I can't reproduce your problem. Please describe your setup (OS,
compiler) and exactly what you did.

The current version of g++ on machine tripped on a couple of
out-of-order declarations in release 1.46. So I fixed these problems
and added your test. Now

git diff r1.46

gives

diff --git a/legacy/C/geodesic.c b/legacy/C/geodesic.c
index 8d9c928..e897e89 100644
--- a/legacy/C/geodesic.c
+++ b/legacy/C/geodesic.c
@@ -187,8 +187,9 @@ static real AngDiff(real x, real y, real* e) {

 static real AngRound(real x) {
   const real z = 1/(real)(16);
+  volatile real y;
   if (x == 0) return 0;
-  volatile real y = fabs(x);
+  y = fabs(x);
   /* The compiler mustn't "simplify" z - (z - y) to y */
   y = y < z ? z - (z - y) : y;
   return x < 0 ? -y : y;
@@ -413,8 +414,8 @@ static void geod_lineinit_int(struct geod_geodesicline* l,
 void geod_lineinit(struct geod_geodesicline* l,
                    const struct geod_geodesic* g,
                    real lat1, real lon1, real azi1, unsigned caps) {
-  azi1 = AngNormalize(azi1);
   real salp1, calp1;
+  azi1 = AngNormalize(azi1);
   /* Guard against underflow in salp0 */
   sincosdx(AngRound(azi1), &salp1, &calp1);
   geod_lineinit_int(l, g, lat1, lon1, azi1, salp1, calp1, caps);
diff --git a/legacy/C/geodtest.c b/legacy/C/geodtest.c
index 990a213..7ec0282 100644
--- a/legacy/C/geodtest.c
+++ b/legacy/C/geodtest.c
@@ -610,6 +610,20 @@ int GeodSolve73() {
   return result;
 }

+int GeodSolve74() {
+  double lat2, lon2, azi2;
+  struct geod_geodesic g;
+  int result = 0;
+  geod_init(&g, wgs84_a, 0);
+  /* geod_init(&g, wgs84_a, wgs84_f); */
+  geod_direct(&g, -3.688856, -104.23828, 0, 4000000,
+              &lat2, &lon2, &azi2);
+  result += assertEquals(lat2, 32.243755, 0.5e-5);
+  result += assertEquals(lon2, -104.23828, 0.5e-5);
+  result += assertEquals(azi2, 0.0, 0.5e-5);
+  return result;
+}
+
 void planimeter(const struct geod_geodesic* g, double points[][2], int N,
                 double* perimeter, double* area) {
   struct geod_polygon p;
@@ -757,6 +771,7 @@ int main() {
   if ((i = GeodSolve67())) {++n; printf("GeodSolve67 fail: %d\n", i);}
   if ((i = GeodSolve71())) {++n; printf("GeodSolve71 fail: %d\n", i);}
   if ((i = GeodSolve73())) {++n; printf("GeodSolve73 fail: %d\n", i);}
+  if ((i = GeodSolve74())) {++n; printf("GeodSolve74 fail: %d\n", i);}
   if ((i = Planimeter0())) {++n; printf("Planimeter0 fail: %d\n", i);}
   if ((i = Planimeter5())) {++n; printf("Planimeter5 fail: %d\n", i);}
   if ((i = Planimeter6())) {++n; printf("Planimeter6 fail: %d\n", i);}

then, on my Linux machine (Fedora 23)

cd legacy/C
mkdir BUILD
cd BUILD
cmake ..
make
make test

gives

Running tests...
Test project /var/tmp/geographiclib/legacy/C/BUILD
    Start 1: geodtest
1/1 Test #1: geodtest .........................   Passed    0.00 sec

100% tests passed, 0 tests failed out of 1

Total Test time (real) =   0.00 sec

[zed]

Sorry, it was my mistake. Now I see that 1.46 works too.