19 if (!(isfinite(_a) && _a > 0))
20 throw GeographicErr(
"Equatorial radius is not positive");
23 throw GeographicErr(
"Gravitational constant is not finite");
27 if (!(isfinite(_omega2) && isfinite(_aomega2)))
28 throw GeographicErr(
"Rotation velocity is not finite");
31 if (!(isfinite(_b) && _b > 0))
32 throw GeographicErr(
"Polar semi-axis is not positive");
35 _ep2 = _e2 / (1 - _e2);
36 real ex2 = _f < 0 ? -_e2 : _ep2;
37 _qQ0 = Qf(ex2, _f < 0);
38 _earth = Geocentric(_a, _f);
39 _eE = _a * sqrt(fabs(_e2));
41 _uU0 = _gGM * atanzz(ex2, _f < 0) / _b + _aomega2 / 3;
42 real P = Hf(ex2, _f < 0) / (6 * _qQ0);
44 _gammae = _gGM / (_a * _b) - (1 + P) * _a * _omega2;
46 _gammap = _gGM / (_a * _a) + 2 * P * _b * _omega2;
49 _k = -_e2 * _gGM / (_a * _b) +
50 _omega2 * (P * (_a + 2 * _b * (1 - _f)) + _a);
52 _fstar = (-_f * _gGM / (_a * _b) + _omega2 * (P * (_a + 2 * _b) + _a)) /
58 Initialize(a, GM, omega, f_J2, geometricp);
83 static const real lg2eps_ = -log2(numeric_limits<real>::epsilon() / 2);
92 int n = x == 0 ? 1 : int(ceil(lg2eps_ / e));
104 return 1/
real(5) + x * atan7series(x);
112 real y = alt ? -x / (1 + x) : x;
113 return !(4 * fabs(y) < 1) ?
114 ((1 + 3/y) * atanzz(x, alt) - 3/y) / (2 * y) :
115 (3 * (3 + y) * atan5series(y) - 1) / 6;
124 real y = alt ? -x / (1 + x) : x;
125 return !(4 * fabs(y) < 1) ?
126 (3 * (1 + 1/y) * (1 - atanzz(x, alt)) - 1) / y :
127 1 - 3 * (1 + y) * atan5series(y);
136 real y = alt ? -x / (1 + x) : x;
137 return !(4 * fabs(y) < 1) ?
138 ((9 + 15/y) * atanzz(x, alt) - 4 - 15/y) / (6 *
Math::sq(y)) :
139 ((25 + 15*y) * atan7series(y) + 3)/10;
148 for (
int j = n; j--;)
151 -3 * e2n * ((1 - n) + 5 * n * _jJ2 / _e2) / ((2 * n + 1) * (2 * n + 3));
161 real& GammaX, real& GammaY, real& GammaZ)
const
166 clam = p != 0 ? X/p : 1,
167 slam = p != 0 ? Y/p : 0,
169 if (_f < 0) swap(p, Z);
176 u = sqrt((Q >= 0 ? (Q + disc) : t2 / (disc - Q)) / 2),
179 sbet = u != 0 ? Z * uE : copysign(sqrt(-Q), Z),
180 cbet = u != 0 ? p * u : p,
181 s = hypot(cbet, sbet);
182 sbet = s != 0 ? sbet/s : 1;
183 cbet = s != 0 ? cbet/s : 0;
187 den = hypot(u, _eE * sbet);
194 bu = _b / (u != 0 || _f < 0 ? u : _eE),
196 q = ((u != 0 || _f < 0 ? Qf(z2, _f < 0) :
Math::pi() / 4) / _qQ0) *
198 qp = _b *
Math::sq(bu) * (u != 0 || _f < 0 ? Hf(z2, _f < 0) : 2) / _qQ0,
199 ang = (
Math::sq(sbet) - 1/real(3)) / 2,
201 Vres = _gGM * (u != 0 || _f < 0 ?
202 atanzz(z2, _f < 0) / u :
203 Math::pi() / (2 * _eE)) + _aomega2 * q * ang,
205 gamu = - (_gGM + (_aomega2 * qp * ang)) * invw /
Math::sq(uE),
206 gamb = _aomega2 * q * sbet * cbet * invw / uE,
208 gamp = t * cbet * gamu - invw * sbet * gamb;
210 GammaX = gamp * clam;
211 GammaY = gamp * slam;
212 GammaZ = invw * sbet * gamu + t * cbet * gamb;
224 real& gammaX, real& gammaY, real& gammaZ)
const {
226 real Ures =
V0(X, Y, Z, gammaX, gammaY, gammaZ) +
Phi(X, Y, fX, fY);
233 real& gammay, real& gammaz)
const {
235 real M[Geocentric::dim2_];
236 _earth.IntForward(lat, 0, h, X, Y, Z, M);
237 real gammaX, gammaY, gammaZ,
238 Ures =
U(X, Y, Z, gammaX, gammaY, gammaZ);
240 gammay = M[1] * gammaX + M[4] * gammaY + M[7] * gammaZ;
241 gammaz = M[2] * gammaX + M[5] * gammaY + M[8] * gammaZ;
246 real omega, real J2) {
250 static const real maxe_ = 1 - numeric_limits<real>::epsilon();
251 static const real eps2_ = sqrt(numeric_limits<real>::epsilon()) / 100;
253 K = 2 *
Math::sq(a * omega) * a / (15 * GM),
255 if (!(GM > 0 && isfinite(K) && K >= 0))
257 if (!(isfinite(J2) && J2 <= J0))
return Math::NaN();
258 if (J2 == J0)
return 1;
266 e2 = fmin(ep2 / (1 + ep2), maxe_);
272 e2a = e2, ep2a = ep2,
275 Q0 = Qf(e2 < 0 ? -e2 : ep2, e2 < 0),
276 h = e2 - f1 * f2 * K / Q0 - 3 * J2,
277 dh = 1 - 3 * f1 * K * QH3f(e2 < 0 ? -e2 : ep2, e2 < 0) /
279 e2 = fmin(e2a - h / dh, maxe_);
280 ep2 = fmax(e2 / (1 - e2), -maxe_);
281 if (fabs(h) < eps2_ || e2 == e2a || ep2 == ep2a)
284 return e2 / (1 + sqrt(1 - e2));
288 real omega, real f) {
290 K = 2 *
Math::sq(a * omega) * a / (15 * GM),
295 return (e2 - K * f1 * f2 / Qf(f < 0 ? -e2 : e2 / f2, f < 0)) / 3;
GeographicLib::Math::real real
#define GEOGRAPHICLIB_PANIC(msg)
Header for GeographicLib::NormalGravity class.
The normal gravity of the earth.
Math::real V0(real X, real Y, real Z, real &GammaX, real &GammaY, real &GammaZ) const
static Math::real FlatteningToJ2(real a, real GM, real omega, real f)
Math::real Phi(real X, real Y, real &fX, real &fY) const
static const NormalGravity & WGS84()
static Math::real J2ToFlattening(real a, real GM, real omega, real J2)
Math::real U(real X, real Y, real Z, real &gammaX, real &gammaY, real &gammaZ) const
Math::real SurfaceGravity(real lat) const
static const NormalGravity & GRS80()
Math::real Gravity(real lat, real h, real &gammay, real &gammaz) const
Namespace for GeographicLib.