12 namespace GeographicLib {
22 , _omega2(
Math::sq(_omega))
23 , _aomega2(
Math::sq(_omega * _a))
28 throw GeographicErr(
"Gravitational constants is not positive");
30 if (_J2 > 0 && Math::isfinite(_J2) && flatp)
32 if (!(_J2 > 0 && Math::isfinite(_J2)) && !flatp)
34 if (!(Math::isfinite(_omega) && _omega != 0))
36 real K = 2 * _aomega2 * _a / (15 * _GM);
39 _ep2 = _e2 / (1 - _e2);
41 _J2 = _e2 * ( 1 - K * sqrt(_e2) / _q0) / 3;
44 for (
int j = 0; j < maxit_; ++j) {
46 real q0 = qf(_e2 / (1 - _e2));
47 _e2 = 3 * _J2 + K * _e2 * sqrt(_e2) / q0;
51 _f = _e2 / (1 + sqrt(1 - _e2));
52 _ep2 = _e2 / (1 - _e2);
58 _U0 = _GM / _E * atan(sqrt(_ep2)) + _aomega2 / 3;
61 _m = _aomega2 * _b / _GM;
63 Q = _m * sqrt(_ep2) * qpf(_ep2) / (3 * _q0),
65 _gammae = _GM / (_a * _b) * G;
66 _gammap = _GM / (_a * _a) * (1 + Q);
68 _k = (_m + 3 * Q / 2 - _e2 * (1 + Q)) / G;
70 _fstar = (_m + 3 * Q / 2 - _f * (1 + Q)) / G;
75 Constants::WGS84_omega<real>(),
76 Constants::WGS84_f<real>(), 0);
80 Constants::GRS80_omega<real>(),
81 0, Constants::GRS80_J2<real>());
92 if (abs(ep2) >
real(0.5))
93 return ((1 + 3 / ep2) * atan(ep) - 3 / ep)/2;
96 for (
int n = 1; ; ++n) {
99 t = (ep2n * n) / ((2 * n + 1) * (2 * n + 3)),
118 if (abs(ep2) >
real(0.5)) {
120 return 3 * (1 + 1 / ep2) * (1 - atan(ep) / ep) - 1;
122 real ep2n = 1, qp = 0;
123 for (
int n = 1; ; ++n) {
126 t = ep2n / ((2 * n + 1) * (2 * n + 3)),
137 Math::real NormalGravity::Jn(
int n)
const throw() {
143 for (
int j = n; j--;)
146 -3 * e2n * (1 - n + 5 * n * _J2 / _e2) / ((2 * n + 1) * (2 * n + 3));
151 phi = lat * Math::degree<real>(),
152 sphi2 = abs(lat) == 90 ? 1 :
Math::sq(sin(phi));
154 return _gammae * (1 + _k * sphi2) / sqrt(1 - _e2 * sphi2);
158 real& GammaX, real& GammaY, real& GammaZ)
171 u = sqrt((Q >= 0 ? (Q + disc) : t2 / (disc - Q)) / 2),
177 cbet = s ? cbet/s : 0;
178 sbet = s ? sbet/s : 1;
186 Vres = (_GM / _E * atan(_E / u)
190 + (_aomega2 * _E * qp
192 gamb = _aomega2 * q * sbet * cbet * invw / uE,
195 GammaX = t * cbet * gamu - invw * sbet * gamb;
196 GammaY = GammaX * slam;
198 GammaZ = invw * sbet * gamu + t * cbet * gamb;
211 real& gammaX, real& gammaY, real& gammaZ)
214 real Ures = V0(X, Y, Z, gammaX, gammaY, gammaZ) + Phi(X, Y, fX, fY);
221 real& gammay, real& gammaz)
224 real M[Geocentric::dim2_];
225 _earth.IntForward(lat, 0, h, X, Y, Z, M);
226 real gammaX, gammaY, gammaZ,
227 Ures = U(X, Y, Z, gammaX, gammaY, gammaZ);
229 gammay = M[1] * gammaX + M[4] * gammaY + M[7] * gammaZ;
230 gammaz = M[2] * gammaX + M[5] * gammaY + M[8] * gammaZ;
Math::real SurfaceGravity(real lat) const
GeographicLib::Math::real real
Math::real Gravity(real lat, real h, real &gammay, real &gammaz) const
The normal gravity of the earth.
static bool isfinite(T x)
static const NormalGravity GRS80
Mathematical functions needed by GeographicLib.
Header for GeographicLib::NormalGravity class.
Exception handling for GeographicLib.
Math::real U(real X, real Y, real Z, real &gammaX, real &gammaY, real &gammaZ) const
static const NormalGravity WGS84
Math::real Phi(real X, real Y, real &fX, real &fY) const
Math::real V0(real X, real Y, real Z, real &GammaX, real &GammaY, real &GammaZ) const