NormalGravity on a sphere, a bug

[Anonymous]

in NormalGravity.cpp, Line 46 ~ 53, an attempt to divide by zero, result is nan

[Charles Karney]

Please try the following patch and let me know if this fixes the problem:

diff --git a/src/NormalGravity.cpp b/src/NormalGravity.cpp
index b673dcd..562feb5 100644
--- a/src/NormalGravity.cpp
+++ b/src/NormalGravity.cpp
@@ -40,20 +40,21 @@ namespace GeographicLib {
         if (!(Math::isfinite(_omega) && _omega != 0))
           throw GeographicErr("Angular velocity is not non-zero");
         if (flatp)
-          _J2 = FlatteningToJ2(a, GM, omega, f);
+          _J2 = FlatteningToJ2(_a, _GM, _omega, _f);
         else
-          _f = J2ToFlattening(a, GM, omega, J2);
+          _f = J2ToFlattening(_a, _GM, _omega, _J2);
       } // else we have a sphere: omega = f = J2 = 0
       _e2 = _f * (2 - _f);
       _ep2 = _e2 / (1 - _e2);
       _q0 = qf(_ep2);
       _earth = Geocentric(_a, _f);
       _b = _a * (1 - _f);
-      _E = a * sqrt(_e2);                               // H+M, Eq 2-54
-      _U0 = _GM / _E * atan(sqrt(_ep2)) + _aomega2 / 3; // H+M, Eq 2-61
+      _E = _a * sqrt(_e2);      // H+M, Eq 2-54
+      // H+M, Eq 2-61
+      _U0 = _GM / (_E ? _E / atan(sqrt(_ep2)) : _b) + _aomega2 / 3;
       // The approximate ratio of the centrifugal acceleration (at the equator)
       // to gravity.
-      _m = _aomega2 * _b / _GM;                         // H+M, Eq 2-70
+      _m = _aomega2 * _b / _GM; // H+M, Eq 2-70
       real
         Q = _m * sqrt(_ep2) * qpf(_ep2) / (3 * _q0),
         G = (1 - _m - Q / 2);

[Anonymous]

that dose not fix the problem, it seems that you still need to fix the function NormalGravity::V0()

  Vres = (_GM / _E * atan(_E / u)
          + _aomega2 * q * (Math::sq(sbet) - 1/real(3)) / 2),

where _E is zero when _f is zero.
the project is really a perfect job; I try to help fix it by self, but I am not familiar with the algorithm, so ... hehe;

[Charles Karney]

Sorry! Here's take 2. Please test this and let me know of any problems.

diff --git a/src/NormalGravity.cpp b/src/NormalGravity.cpp
index b673dcd..f924d86 100644
--- a/src/NormalGravity.cpp
+++ b/src/NormalGravity.cpp
@@ -30,7 +30,7 @@ namespace GeographicLib {
       if (!(Math::isfinite(_a) && _a > 0))
         throw GeographicErr("Major radius is not positive");
       if (!(Math::isfinite(_GM) && _GM > 0))
-        throw GeographicErr("Gravitational constants is not positive");
+        throw GeographicErr("Gravitational constant is not positive");
       if (!(_omega == 0 && _f == 0 && _J2 == 0)) {
         bool flatp = _f > 0 && Math::isfinite(_f);
         if (_J2 > 0 && Math::isfinite(_J2) && flatp)
@@ -40,22 +40,23 @@ namespace GeographicLib {
         if (!(Math::isfinite(_omega) && _omega != 0))
           throw GeographicErr("Angular velocity is not non-zero");
         if (flatp)
-          _J2 = FlatteningToJ2(a, GM, omega, f);
+          _J2 = FlatteningToJ2(_a, _GM, _omega, _f);
         else
-          _f = J2ToFlattening(a, GM, omega, J2);
+          _f = J2ToFlattening(_a, _GM, _omega, _J2);
       } // else we have a sphere: omega = f = J2 = 0
       _e2 = _f * (2 - _f);
       _ep2 = _e2 / (1 - _e2);
       _q0 = qf(_ep2);
       _earth = Geocentric(_a, _f);
       _b = _a * (1 - _f);
-      _E = a * sqrt(_e2);                               // H+M, Eq 2-54
-      _U0 = _GM / _E * atan(sqrt(_ep2)) + _aomega2 / 3; // H+M, Eq 2-61
+      _E = _a * sqrt(_e2);      // H+M, Eq 2-54
+      // H+M, Eq 2-61
+      _U0 = _GM / (_E ? _E / atan(sqrt(_ep2)) : _b) + _aomega2 / 3;
       // The approximate ratio of the centrifugal acceleration (at the equator)
       // to gravity.
-      _m = _aomega2 * _b / _GM;                         // H+M, Eq 2-70
+      _m = _aomega2 * _b / _GM; // H+M, Eq 2-70
       real
-        Q = _m * sqrt(_ep2) * qpf(_ep2) / (3 * _q0),
+        Q = _m * (_q0 ? sqrt(_ep2) * qpf(_ep2) / (3 * _q0) : 1),
         G = (1 - _m - Q / 2);
       _gammae = _GM / (_a * _b) * G;       // H+M, Eq 2-73
       _gammap = _GM / (_a * _a) * (1 + Q); // H+M, Eq 2-74
@@ -179,14 +180,14 @@ namespace GeographicLib {
       invw = uE / Math::hypot(u, _E * sbet), // H+M, Eq 2-63
       ep = _E/u,
       ep2 = Math::sq(ep),
-      q = qf(ep2) / _q0,
+      q = _q0 ? qf(ep2) / _q0 : pow(_a / u, 3),
       qp = qpf(ep2) / _q0,
       // H+M, Eqs 2-62 + 6-9, but omitting last (rotational) term .
-      Vres = (_GM / _E * atan(_E / u)
+      Vres = (_GM / (_E ? _E / atan(_E / u) : u)
               + _aomega2 * q * (Math::sq(sbet) - 1/real(3)) / 2),
       // H+M, Eq 6-10
       gamu = - invw * (_GM
-                       + (_aomega2 * _E * qp
+                       + (_aomega2 * (_q0 ? _E * qp : 3 * q * u)
                           * (Math::sq(sbet) - 1/real(3)) / 2)) / Math::sq(uE),
       gamb = _aomega2 * q * sbet * cbet * invw / uE,
       t = u * invw / uE;

[Htallon]

I think the patch should fix the problem: when I use this:

NormalGravity grav(6371009, Constants::WGS84_GM(),
            0, 0, 0);
double gammay, gammaz;
grav.Gravity(30, 0, gammay, gammaz);

the result gammaz is -9.82022 at different latitudes, with the height=0; it is reasonable.

I have another question, maybe not necessary: why can't I compute NormalGravity on the nonrotating ellipsoid or rotating sphere? that is when I use in such a way:

NormalGravity grav(Constants::WGS84_a(), Constants::WGS84_GM(),
                0, Constants::WGS84_f(), 0);

and

NormalGravity grav(6371009, Constants::WGS84_GM(),
            Constants::WGS84_omega(), 0, 0);

they both throw error or exception. So I wonder if I need to use other methods of NormalGravity to do such computation or I don't need to do such computation in fact, maybe the result doesn't matter.

If I want to model the gravity on the nonrotating ellipsoid and rotating sphere each, how to use this lib in an uniform way?

Hope you can understand what I say for bearing my poor English.
Thanks

[Charles Karney]

The normal gravity is a vacuum gravitational potential V0 such that the
sum of V0 and the centrifugal potential Phi is a constant on an
ellipsoid of revolution. The solution to this problem, due to
Somigliana and described in Heiskanen and Moritz (1967), is such that as
Omega approaches 0, the flattening, f, must also approach 0 (as
Omega^2). It may be that there's some way of taking the limits to
remove this restriction (and also to allow for prolate ellipsoids, with
f < 0). However, the NormalGravity class does not allow this at
present.

Now to answer your specific request for the gravity of a rotating sphere
(assumed to be of constant density). This is simple, just add the 1/r^2
gravity and the centrifugal force. A similar calculation can be done
for a non-rotating ellipsoid of constant density (I believe Laplace
computed the gravity in this case). However neither of these is a
normal gravity solution, because the total potential is not constant
on the surface in question.

Incidentally, in the normal gravity problem, nothing is said about the
mass distribution within the ellipsoid. (I'm pretty sure that a
constant density doesn't yield the desired potential.)