14 # pragma warning (disable: 4127) 27 , _omega2(
Math::sq(_omega))
28 , _aomega2(
Math::sq(_omega * _a))
33 throw GeographicErr(
"Gravitational constant is not positive");
34 if (!(_omega == 0 && _f == 0 && _J2 == 0)) {
36 if (_J2 > 0 && Math::isfinite(_J2) && flatp)
38 if (!(_J2 > 0 && Math::isfinite(_J2)) && !flatp)
40 if (!(Math::isfinite(_omega) && _omega != 0))
48 _ep2 = _e2 / (1 - _e2);
54 _U0 = _GM / (_E ? _E / atan(sqrt(_ep2)) : _b) + _aomega2 / 3;
57 _m = _aomega2 * _b / _GM;
59 Q = _m * (_q0 ? sqrt(_ep2) * qpf(_ep2) / (3 * _q0) : 1),
61 _gammae = _GM / (_a * _b) * G;
62 _gammap = _GM / (_a * _a) * (1 + Q);
64 _k = (_m + 3 * Q / 2 - _e2 * (1 + Q)) / G;
66 _fstar = (_m + 3 * Q / 2 - _f * (1 + Q)) / G;
87 if (abs(x) >= real(0.5)) {
88 real y = sqrt(abs(x)), x2 =
Math::sq(x);
89 return ((x > 0 ? atan(y) :
Math::atanh(y))- y * (1 - x / 3 + x2 / 5)) /
93 for (
int n = 7; ; n += 2) {
106 {
return 1/real(5) + x * atan7(x); }
116 return sqrt(ep2) * ep2 * (3 * (3 + ep2) * atan5(ep2) - 1) / 6;
122 return sqrt(ep2) * (5 - 3 * (1 + ep2) * (1 + 5 * ep2 * atan7(ep2))) /
134 return ep2 * (1 - 3 * (1 + ep2) * atan5(ep2));
143 for (
int j = n; j--;)
146 -3 * e2n * ((1 - n) + 5 * n * _J2 / _e2) / ((2 * n + 1) * (2 * n + 3));
152 return _gammae * (1 + _k *
Math::sq(sphi)) / sqrt(1 - _e2 *
Math::sq(sphi));
156 real& GammaX, real& GammaY, real& GammaZ)
169 u = sqrt((Q >= 0 ? (Q + disc) : t2 / (disc - Q)) / 2),
175 cbet = s ? cbet/s : 0;
176 sbet = s ? sbet/s : 1;
181 q = _q0 ? qf(ep2) / _q0 : pow(_a / u, 3),
184 Vres = (_GM / (_E ? _E / atan(_E / u) : u)
185 + _aomega2 * q * (
Math::sq(sbet) - 1/real(3)) / 2),
188 + (_aomega2 * (_q0 ? _E * qp : 3 * q * u)
190 gamb = _aomega2 * q * sbet * cbet * invw / uE,
193 GammaX = t * cbet * gamu - invw * sbet * gamb;
194 GammaY = GammaX * slam;
196 GammaZ = invw * sbet * gamu + t * cbet * gamb;
209 real& gammaX, real& gammaY, real& gammaZ)
212 real Ures =
V0(X, Y, Z, gammaX, gammaY, gammaZ) +
Phi(X, Y, fX, fY);
219 real& gammay, real& gammaz)
222 real M[Geocentric::dim2_];
223 _earth.IntForward(lat, 0, h, X, Y, Z, M);
224 real gammaX, gammaY, gammaZ,
225 Ures =
U(X, Y, Z, gammaX, gammaY, gammaZ);
227 gammay = M[1] * gammaX + M[4] * gammaY + M[7] * gammaZ;
228 gammaz = M[2] * gammaX + M[5] * gammaY + M[8] * gammaZ;
233 real omega, real J2) {
235 K = 2 *
Math::sq(a * omega) * a / (15 * GM),
243 h = e2 * (1 - sqrt(e2) * K / q0) - 3 * J2,
244 dh = 1 - sqrt(e2) * K * (3 * q0 - 2 * e2 * dq0) / (2 *
Math::sq(q0)),
250 return e2 / (1 + sqrt(1 - e2));
254 real omega, real f) {
256 K = 2 *
Math::sq(a * omega) * a / (15 * GM),
258 q0 = qf(e2 / (1 - e2));
259 return e2 * (1 - K * sqrt(e2) / q0) / 3;
Math::real SurfaceGravity(real lat) const
Math::real Gravity(real lat, real h, real &gammay, real &gammaz) const
The normal gravity of the earth.
static bool isfinite(T x)
Mathematical functions needed by GeographicLib.
static Math::real FlatteningToJ2(real a, real GM, real omega, real f)
Namespace for GeographicLib.
static Math::real J2ToFlattening(real a, real GM, real omega, real J2)
static const NormalGravity & GRS80()
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
#define GEOGRAPHICLIB_PANIC
Math::real V0(real X, real Y, real Z, real &GammaX, real &GammaY, real &GammaZ) const