GeographicLib: NormalGravity.cpp Source File (original) (raw)

14 using namespace std;

16 void NormalGravity::Initialize(real a, real GM, real omega, real f_J2,

17 bool geometricp) {

18 _a = a;

19 if (!(isfinite(_a) && _a > 0))

20 throw GeographicErr("Equatorial radius is not positive");

21 _gGM = GM;

22 if (!isfinite(_gGM))

23 throw GeographicErr("Gravitational constant is not finite");

24 _omega = omega;

26 _aomega2 = Math::sq(_omega * _a);

27 if (!(isfinite(_omega2) && isfinite(_aomega2)))

28 throw GeographicErr("Rotation velocity is not finite");

29 _f = geometricp ? f_J2 : J2ToFlattening(_a, _gGM, _omega, f_J2);

30 _b = _a * (1 - _f);

31 if (!(isfinite(_b) && _b > 0))

32 throw GeographicErr("Polar semi-axis is not positive");

33 _jJ2 = geometricp ? FlatteningToJ2(_a, _gGM, _omega, f_J2) : f_J2;

34 _e2 = _f * (2 - _f);

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)) /

53 _gammae;

54 }

57 bool geometricp) {

58 Initialize(a, GM, omega, f_J2, geometricp);

59 }

83 static const real lg2eps_ = -log2(numeric_limits::epsilon() / 2);

84 int e;

85 (void) frexp(x, &e);

86 e = max(-e, 1);

92 int n = x == 0 ? 1 : int(ceil(lg2eps_ / e));

94 while (n--)

96 return v;

97 }

100

101

102

103

104 return 1/real(5) + x * atan7series(x);

105 }

106

108

109

110

111

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;

116 }

117

119

120

121

122

123

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);

128 }

129

131

132

133

134

135

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;

140 }

141

142 Math::real NormalGravity::Jn(int n) const {

143

144 if (n & 1 || n < 0)

145 return 0;

146 n /= 2;

147 real e2n = 1;

148 for (int j = n; j--;)

149 e2n *= -_e2;

150 return

151 -3 * e2n * ((1 - n) + 5 * n * _jJ2 / _e2) / ((2 * n + 1) * (2 * n + 3));

152 }

153

156

157 return (_gammae + _k * Math::sq(sphi)) / sqrt(1 - _e2 * Math::sq(sphi));

158 }

159

161 real& GammaX, real& GammaY, real& GammaZ) const

162 {

163

164 real

165 p = hypot(X, Y),

166 clam = p != 0 ? X/p : 1,

167 slam = p != 0 ? Y/p : 0,

168 r = hypot(p, Z);

169 if (_f < 0) swap(p, Z);

170 real

173 disc = sqrt(Math::sq(Q) + t2),

174

175

176 u = sqrt((Q >= 0 ? (Q + disc) : t2 / (disc - Q)) / 2),

177 uE = hypot(u, _eE),

178

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;

184 real

185 z = _eE/u,

187 den = hypot(u, _eE * sbet);

188 if (_f < 0) {

189 swap(sbet, cbet);

191 }

192 real

193 invw = uE / den,

194 bu = _b / (u != 0 || _f < 0 ? u : _eE),

195

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,

200

201 Vres = _gGM * (u != 0 || _f < 0 ?

202 atanzz(z2, _f < 0) / u :

203 Math::pi() / (2 * _eE)) + _aomega2 * q * ang,

204

205 gamu = - (_gGM + (_aomega2 * qp * ang)) * invw / Math::sq(uE),

206 gamb = _aomega2 * q * sbet * cbet * invw / uE,

207 t = u * invw / uE,

208 gamp = t * cbet * gamu - invw * sbet * gamb;

209

210 GammaX = gamp * clam;

211 GammaY = gamp * slam;

212 GammaZ = invw * sbet * gamu + t * cbet * gamb;

213 return Vres;

214 }

215

217 fX = _omega2 * X;

218 fY = _omega2 * Y;

219

221 }

222

224 real& gammaX, real& gammaY, real& gammaZ) const {

225 real fX, fY;

226 real Ures = V0(X, Y, Z, gammaX, gammaY, gammaZ) + Phi(X, Y, fX, fY);

227 gammaX += fX;

228 gammaY += fY;

229 return Ures;

230 }

231

233 real& gammay, real& gammaz) const {

234 real X, Y, Z;

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);

239

240 gammay = M[1] * gammaX + M[4] * gammaY + M[7] * gammaZ;

241 gammaz = M[2] * gammaX + M[5] * gammaY + M[8] * gammaZ;

242 return Ures;

243 }

244

246 real omega, real J2) {

247

248

249

250 static const real maxe_ = 1 - numeric_limits::epsilon();

251 static const real eps2_ = sqrt(numeric_limits::epsilon()) / 100;

252 real

253 K = 2 * Math::sq(a * omega) * a / (15 * GM),

254 J0 = (1 - 4 * K / Math::pi()) / 3;

255 if (!(GM > 0 && isfinite(K) && K >= 0))

257 if (!(isfinite(J2) && J2 <= J0)) return Math::NaN();

258 if (J2 == J0) return 1;

259

260

261

262

263 real

265 -maxe_),

266 e2 = fmin(ep2 / (1 + ep2), maxe_);

267 for (int j = 0;

268 j < maxit_ ||

270 ++j) {

271 real

272 e2a = e2, ep2a = ep2,

273 f2 = 1 - e2,

274 f1 = sqrt(f2),

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)

282 break;

283 }

284 return e2 / (1 + sqrt(1 - e2));

285 }

286

288 real omega, real f) {

289 real

290 K = 2 * Math::sq(a * omega) * a / (15 * GM),

291 f1 = 1 - f,

293 e2 = f * (2 - f);

294

295 return (e2 - K * f1 * f2 / Qf(f < 0 ? -e2 : e2 / f2, f < 0)) / 3;

296 }

297

298}

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

Definition NormalGravity.cpp:160

static Math::real FlatteningToJ2(real a, real GM, real omega, real f)

Definition NormalGravity.cpp:287

Math::real Phi(real X, real Y, real &fX, real &fY) const

Definition NormalGravity.cpp:216

static const NormalGravity & WGS84()

Definition NormalGravity.cpp:61

static Math::real J2ToFlattening(real a, real GM, real omega, real J2)

Definition NormalGravity.cpp:245

Math::real U(real X, real Y, real Z, real &gammaX, real &gammaY, real &gammaZ) const

Definition NormalGravity.cpp:223

Math::real SurfaceGravity(real lat) const

Definition NormalGravity.cpp:154

static const NormalGravity & GRS80()

Definition NormalGravity.cpp:69

Math::real Gravity(real lat, real h, real &gammay, real &gammaz) const

Definition NormalGravity.cpp:232

Namespace for GeographicLib.

void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)