14 # pragma warning (disable: 4127) 25 throw GeographicErr(
"Equatorial radius is not positive");
28 throw GeographicErr(
"Gravitational constant is not finite");
33 throw GeographicErr(
"Rotation velocity is not finite");
34 _f = geometricp ? f_J2 : J2ToFlattening(_a, _GM, _omega, f_J2);
37 throw GeographicErr(
"Polar semi-axis is not positive");
38 _J2 = geometricp ? FlatteningToJ2(_a, _GM, _omega, f_J2) : f_J2;
40 _ep2 = _e2 / (1 - _e2);
41 real ex2 = _f < 0 ? -_e2 : _ep2;
42 _Q0 = Qf(ex2, _f < 0);
43 _earth = Geocentric(_a, _f);
44 _E = _a * sqrt(abs(_e2));
46 _U0 = _GM * atanzz(ex2, _f < 0) / _b + _aomega2 / 3;
47 real P = Hf(ex2, _f < 0) / (6 * _Q0);
49 _gammae = _GM / (_a * _b) - (1 + P) * _a * _omega2;
51 _gammap = _GM / (_a * _a) + 2 * P * _b * _omega2;
54 _k = -_e2 * _GM / (_a * _b) +
55 _omega2 * (P * (_a + 2 * _b * (1 - _f)) + _a);
57 _fstar = (-_f * _GM / (_a * _b) + _omega2 * (P * (_a + 2 * _b) + _a)) /
63 Initialize(a, GM, omega, f_J2, geometricp);
68 throw GeographicErr(
"Gravitational constant is not positive");
70 if (!(omega == 0 && f == 0 && J2 == 0)) {
80 Initialize(a, GM, omega, geometricp ? f : J2, geometricp);
105 static const real lg2eps_ =
106 -log(numeric_limits<real>::epsilon() / 2) / log(
real(2));
115 int n = int(ceil(lg2eps_ / e));
127 return 1/
real(5) + x * atan7series(x);
135 real y = alt ? -x / (1 + x) : x;
136 return !(4 * abs(y) < 1) ?
137 ((1 + 3/y) * atanzz(x, alt) - 3/y) / (2 * y) :
138 (3 * (3 + y) * atan5series(y) - 1) / 6;
147 real y = alt ? -x / (1 + x) : x;
148 return !(4 * abs(y) < 1) ?
149 (3 * (1 + 1/y) * (1 - atanzz(x, alt)) - 1) / y :
150 1 - 3 * (1 + y) * atan5series(y);
159 real y = alt ? -x / (1 + x) : x;
160 return !(4 * abs(y) < 1) ?
161 ((9 + 15/y) * atanzz(x, alt) - 4 - 15/y) / (6 *
Math::sq(y)) :
162 ((25 + 15*y) * atan7series(y) + 3)/10;
171 for (
int j = n; j--;)
174 -3 * e2n * ((1 - n) + 5 * n * _J2 / _e2) / ((2 * n + 1) * (2 * n + 3));
184 real& GammaX, real& GammaY, real& GammaZ)
const 189 clam = p != 0 ? X/p : 1,
190 slam = p != 0 ? Y/p : 0,
192 if (_f < 0)
swap(p, Z);
199 u = sqrt((Q >= 0 ? (Q + disc) : t2 / (disc - Q)) / 2),
203 cbet = u != 0 ? p * u : p,
205 sbet = s != 0 ? sbet/s : 1;
206 cbet = s != 0 ? cbet/s : 0;
217 bu = _b / (u != 0 || _f < 0 ? u : _E),
219 q = ((u != 0 || _f < 0 ? Qf(z2, _f < 0) :
Math::pi() / 4) / _Q0) *
221 qp = _b *
Math::sq(bu) * (u != 0 || _f < 0 ? Hf(z2, _f < 0) : 2) / _Q0,
224 Vres = _GM * (u != 0 || _f < 0 ?
225 atanzz(z2, _f < 0) / u :
226 Math::pi() / (2 * _E)) + _aomega2 * q * ang,
228 gamu = - (_GM + (_aomega2 * qp * ang)) * invw /
Math::sq(uE),
229 gamb = _aomega2 * q * sbet * cbet * invw / uE,
231 gamp = t * cbet * gamu - invw * sbet * gamb;
233 GammaX = gamp * clam;
234 GammaY = gamp * slam;
235 GammaZ = invw * sbet * gamu + t * cbet * gamb;
247 real& gammaX, real& gammaY, real& gammaZ)
const {
249 real Ures = V0(X, Y, Z, gammaX, gammaY, gammaZ) + Phi(X, Y, fX, fY);
256 real& gammay, real& gammaz)
const {
258 real M[Geocentric::dim2_];
259 _earth.IntForward(lat, 0, h, X, Y, Z, M);
260 real gammaX, gammaY, gammaZ,
261 Ures = U(X, Y, Z, gammaX, gammaY, gammaZ);
263 gammay = M[1] * gammaX + M[4] * gammaY + M[7] * gammaZ;
264 gammaz = M[2] * gammaX + M[5] * gammaY + M[8] * gammaZ;
269 real omega, real J2) {
273 static const real maxe_ = 1 - numeric_limits<real>::epsilon();
274 static const real eps2_ = sqrt(numeric_limits<real>::epsilon()) / 100;
276 K = 2 *
Math::sq(a * omega) * a / (15 * GM),
281 if (J2 == J0)
return 1;
289 e2 = min(ep2 / (1 + ep2), maxe_);
292 e2a = e2, ep2a = ep2,
295 Q0 = Qf(e2 < 0 ? -e2 : ep2, e2 < 0),
296 h = e2 - f1 * f2 * K / Q0 - 3 * J2,
297 dh = 1 - 3 * f1 * K * QH3f(e2 < 0 ? -e2 : ep2, e2 < 0) /
299 e2 = min(e2a - h / dh, maxe_);
300 ep2 = max(e2 / (1 - e2), -maxe_);
301 if (abs(h) < eps2_ || e2 == e2a || ep2 == ep2a)
304 return e2 / (1 + sqrt(1 - e2));
308 real omega, real f) {
310 K = 2 *
Math::sq(a * omega) * a / (15 * GM),
315 return (e2 - K * f1 * f2 / Qf(f < 0 ? -e2 : e2 / f2, f < 0)) / 3;
GeographicLib::Math::real real
The normal gravity of the earth.
static bool isfinite(T x)
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()
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
Math::real SurfaceGravity(real lat) const
static T copysign(T x, T y)
Header for GeographicLib::NormalGravity class.
Exception handling for GeographicLib.
Math::real Gravity(real lat, real h, real &gammay, real &gammaz) const
static const NormalGravity & WGS84()
Math::real V0(real X, real Y, real Z, real &GammaX, real &GammaY, real &GammaZ) const
Math::real Phi(real X, real Y, real &fX, real &fY) const
Math::real U(real X, real Y, real Z, real &gammaX, real &gammaY, real &gammaZ) const
#define GEOGRAPHICLIB_PANIC