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

11#include

12#include

17#if !defined(GEOGRAPHICLIB_DATA)

18# if defined(_WIN32)

19# define GEOGRAPHICLIB_DATA "C:/ProgramData/GeographicLib"

20# else

21# define GEOGRAPHICLIB_DATA "/usr/local/share/GeographicLib"

22# endif

23#endif

25#if !defined(GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME)

26# define GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME "egm96"

27#endif

29#if defined(_MSC_VER)

31# pragma warning (disable: 4996)

32#endif

36 using namespace std;

38 GravityModel::GravityModel(const std::string& name, const std::string& path,

39 int Nmax, int Mmax)

40 : _name(name)

41 , _dir(path)

42 , _description("NONE")

43 , _date("UNKNOWN")

44 , _amodel(Math::NaN())

45 , _gGMmodel(Math::NaN())

46 , _zeta0(0)

47 , _corrmult(1)

48 , _nmx(-1)

49 , _mmx(-1)

51 {

52 if (_dir.empty())

54 bool truncate = Nmax >= 0 || Mmax >= 0;

55 if (truncate) {

56 if (Nmax >= 0 && Mmax < 0) Mmax = Nmax;

57 if (Nmax < 0) Nmax = numeric_limits::max();

58 if (Mmax < 0) Mmax = numeric_limits::max();

59 }

60 ReadMetadata(_name);

61 {

62 string coeff = _filename + ".cof";

63 ifstream coeffstr(coeff.c_str(), ios::binary);

64 if (!coeffstr.good())

66 char id[idlength_ + 1];

67 coeffstr.read(id, idlength_);

68 if (!coeffstr.good())

70 id[idlength_] = '\0';

71 if (_id != string(id))

72 throw GeographicErr("ID mismatch: " + _id + " vs " + id);

73 int N, M;

74 if (truncate) { N = Nmax; M = Mmax; }

76 if (!(N >= 0 && M >= 0))

77 throw GeographicErr("Degree and order must be at least 0");

78 if (_cCx[0] != 0)

79 throw GeographicErr("The degree 0 term should be zero");

80 _cCx[0] = 1;

81 _gravitational = SphericalHarmonic(_cCx, _sSx, N, N, M, _amodel, _norm);

82 if (truncate) { N = Nmax; M = Mmax; }

84 if (N < 0) {

85 N = M = 0;

86 _cCC.resize(1, real(0));

87 }

88 _cCC[0] += _zeta0 / _corrmult;

89 _correction = SphericalHarmonic(_cCC, _cCS, N, N, M, real(1), _norm);

90 int pos = int(coeffstr.tellg());

91 coeffstr.seekg(0, ios::end);

92 if (pos != coeffstr.tellg())

94 }

100 real mult = _earth._gGM / _gGMmodel;

101 real amult = Math::sq(_earth._a / _amodel);

102

103

104

105 _zonal.clear(); _zonal.push_back(1);

106 _dzonal0 = (_earth.MassConstant() - _gGMmodel) / _gGMmodel;

107 for (int n = 2; n <= nmx; n += 2) {

108

109

110

111

112

113 mult *= amult;

114 real

115 r = _cCx[n],

116 s = - mult * _earth.Jn(n) / sqrt(real(2 * n + 1)),

117 t = r - s;

118 if (t == r)

119 break;

120 _zonal.push_back(0);

121 _zonal.push_back(s);

122 }

123 int nmx1 = int(_zonal.size()) - 1;

127 _zonal,

128 _zonal,

129 nmx1, nmx1, 0,

130 _amodel,

132 }

133

134 void GravityModel::ReadMetadata(const string& name) {

135 const char* spaces = " \t\n\v\f\r";

136 _filename = _dir + "/" + name + ".egm";

137 ifstream metastr(_filename.c_str());

138 if (!metastr.good())

140 string line;

141 getline(metastr, line);

142 if (!(line.size() >= 6 && line.substr(0,5) == "EGMF-"))

143 throw GeographicErr(_filename + " does not contain EGMF-n signature");

144 string::size_type n = line.find_first_of(spaces, 5);

145 if (n != string::npos)

146 n -= 5;

147 string version(line, 5, n);

148 if (version != "1")

149 throw GeographicErr("Unknown version in " + _filename + ": " + version);

150 string key, val;

151 real a = Math::NaN(), GM = a, omega = a, f = a, J2 = a;

152 while (getline(metastr, line)) {

154 continue;

155

156 if (key == "Name")

157 _name = val;

158 else if (key == "Description")

159 _description = val;

160 else if (key == "ReleaseDate")

161 _date = val;

162 else if (key == "ModelRadius")

163 _amodel = Utility::val(val);

164 else if (key == "ModelMass")

165 _gGMmodel = Utility::val(val);

166 else if (key == "AngularVelocity")

167 omega = Utility::val(val);

168 else if (key == "ReferenceRadius")

169 a = Utility::val(val);

170 else if (key == "ReferenceMass")

171 GM = Utility::val(val);

172 else if (key == "Flattening")

173 f = Utility::fract(val);

174 else if (key == "DynamicalFormFactor")

175 J2 = Utility::fract(val);

176 else if (key == "HeightOffset")

177 _zeta0 = Utility::fract(val);

178 else if (key == "CorrectionMultiplier")

179 _corrmult = Utility::fract(val);

180 else if (key == "Normalization") {

181 if (val == "FULL" || val == "Full" || val == "full")

183 else if (val == "SCHMIDT" || val == "Schmidt" || val == "schmidt")

185 else

186 throw GeographicErr("Unknown normalization " + val);

187 } else if (key == "ByteOrder") {

188 if (val == "Big" || val == "big")

189 throw GeographicErr("Only little-endian ordering is supported");

190 else if (!(val == "Little" || val == "little"))

191 throw GeographicErr("Unknown byte ordering " + val);

192 } else if (key == "ID")

193 _id = val;

194

195 }

196

197 if (!(isfinite(_amodel) && _amodel > 0))

198 throw GeographicErr("Model radius must be positive");

199 if (!(isfinite(_gGMmodel) && _gGMmodel > 0))

200 throw GeographicErr("Model mass constant must be positive");

201 if (!(isfinite(_corrmult) && _corrmult > 0))

202 throw GeographicErr("Correction multiplier must be positive");

203 if (!(isfinite(_zeta0)))

204 throw GeographicErr("Height offset must be finite");

205 if (int(_id.size()) != idlength_)

206 throw GeographicErr("Invalid ID");

207 if (isfinite(f) && isfinite(J2))

208 throw GeographicErr("Cannot specify both f and J2");

209 _earth = NormalGravity(a, GM, omega,

210 isfinite(f) ? f : J2, isfinite(f));

211 }

212

215 bool gradp, bool correct) const {

216

217

218

219 if (_dzonal0 == 0)

220

221 correct = false;

222 real T, invR = correct ? 1 / hypot(hypot(X, Y), Z) : 1;

223 if (gradp) {

224

225 deltaX = deltaY = deltaZ = 0;

226 T = _disturbing(-1, X, Y, Z, deltaX, deltaY, deltaZ);

227 real f = _gGMmodel / _amodel;

228 deltaX *= f;

229 deltaY *= f;

230 deltaZ *= f;

231 if (correct) {

232 invR = _gGMmodel * _dzonal0 * invR * invR * invR;

233 deltaX += X * invR;

234 deltaY += Y * invR;

235 deltaZ += Z * invR;

236 }

237 } else

238 T = _disturbing(-1, X, Y, Z);

239 T = (T / _amodel - (correct ? _dzonal0 : 0) * invR) * _gGMmodel;

240 return T;

241 }

242

244 real& GX, real& GY, real& GZ) const {

245 real

246 Vres = _gravitational(X, Y, Z, GX, GY, GZ),

247 f = _gGMmodel / _amodel;

248 Vres *= f;

249 GX *= f;

250 GY *= f;

251 GZ *= f;

252 return Vres;

253 }

254

256 real& gX, real& gY, real& gZ) const {

257 real fX, fY,

258 Wres = V(X, Y, Z, gX, gY, gZ) + _earth.Phi(X, Y, fX, fY);

259 gX += fX;

260 gY += fY;

261 return Wres;

262 }

263

265 real& Dg01, real& xi, real& eta) const {

266 real X, Y, Z, M[Geocentric::dim2_];

267 _earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);

268 real

269 deltax, deltay, deltaz,

270 T = InternalT(X, Y, Z, deltax, deltay, deltaz, true, false),

271 clam = M[3], slam = -M[0],

272 P = hypot(X, Y),

273 R = hypot(P, Z),

274

275 ctheta = R != 0 ? P / R : M[7],

276 stheta = R != 0 ? Z / R : M[8];

277

278 real MC[Geocentric::dim2_];

279 Geocentric::Rotation(stheta, ctheta, slam, clam, MC);

280 Geocentric::Unrotate(MC, deltax, deltay, deltaz, deltax, deltay, deltaz);

281

282 Dg01 = - deltaz - 2 * T / R;

283 real gammaX, gammaY, gammaZ;

284 _earth.U(X, Y, Z, gammaX, gammaY, gammaZ);

285 real gamma = hypot( hypot(gammaX, gammaY), gammaZ);

288 }

289

291 {

292 real X, Y, Z;

293 _earth.Earth().IntForward(lat, lon, 0, X, Y, Z, NULL);

294 real

296 dummy,

297 T = InternalT(X, Y, Z, dummy, dummy, dummy, false, false),

298 invR = 1 / hypot(hypot(X, Y), Z),

299 correction = _corrmult * _correction(invR * X, invR * Y, invR * Z);

300

301 return T/gamma0 + correction;

302 }

303

305 real& gx, real& gy, real& gz) const {

306 real X, Y, Z, M[Geocentric::dim2_];

307 _earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);

308 real Wres = W(X, Y, Z, gx, gy, gz);

309 Geocentric::Unrotate(M, gx, gy, gz, gx, gy, gz);

310 return Wres;

311 }

313 real& deltax, real& deltay,

314 real& deltaz) const {

315 real X, Y, Z, M[Geocentric::dim2_];

316 _earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);

317 real Tres = InternalT(X, Y, Z, deltax, deltay, deltaz, true, true);

318 Geocentric::Unrotate(M, deltax, deltay, deltaz, deltax, deltay, deltaz);

319 return Tres;

320 }

321

323 if (h != 0)

324

325 caps &= ~(CAP_GAMMA0 | CAP_C);

326 real X, Y, Z, M[Geocentric::dim2_];

327 _earth.Earth().IntForward(lat, 0, h, X, Y, Z, M);

328

329 real

330 invR = 1 / hypot(X, Z),

331 gamma0 = (caps & CAP_GAMMA0 ?_earth.SurfaceGravity(lat)

333 fx, fy, fz, gamma;

334 if (caps & CAP_GAMMA) {

335 _earth.U(X, Y, Z, fx, fy, fz);

336 gamma = hypot(fx, fz);

337 } else

339 _earth.Phi(X, Y, fx, fy);

341 _earth._a, _earth._f, lat, h, Z, X, M[7], M[8],

342 _amodel, _gGMmodel, _dzonal0, _corrmult,

343 gamma0, gamma, fx,

344 caps & CAP_G ?

345 _gravitational.Circle(X, Z, true) :

347

348 caps & CAP_T ?

349 _disturbing.Circle(-1, X, Z, (caps&CAP_DELTA) != 0) :

351 caps & CAP_C ?

352 _correction.Circle(invR * X, invR * Z, false) :

354 }

355

357 string path;

358 char* gravitypath = getenv("GEOGRAPHICLIB_GRAVITY_PATH");

359 if (gravitypath)

360 path = string(gravitypath);

361 if (!path.empty())

362 return path;

363 char* datapath = getenv("GEOGRAPHICLIB_DATA");

364 if (datapath)

365 path = string(datapath);

366 return (!path.empty() ? path : string(GEOGRAPHICLIB_DATA)) + "/gravity";

367 }

368

370 string name;

371 char* gravityname = getenv("GEOGRAPHICLIB_GRAVITY_NAME");

372 if (gravityname)

373 name = string(gravityname);

375 }

376

377}

GeographicLib::Math::real real

#define GEOGRAPHICLIB_DATA

Header for GeographicLib::GravityCircle class.

#define GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME

Definition GravityModel.cpp:26

Header for GeographicLib::GravityModel class.

Header for GeographicLib::SphericalEngine class.

Header for GeographicLib::Utility class.

Spherical harmonic sums for a circle.

Exception handling for GeographicLib.

Gravity on a circle of latitude.

friend class GravityCircle

Math::real Disturbance(real lat, real lon, real h, real &deltax, real &deltay, real &deltaz) const

Definition GravityModel.cpp:312

static std::string DefaultGravityPath()

Definition GravityModel.cpp:356

void SphericalAnomaly(real lat, real lon, real h, real &Dg01, real &xi, real &eta) const

Definition GravityModel.cpp:264

GravityCircle Circle(real lat, real h, unsigned caps=ALL) const

Definition GravityModel.cpp:322

static std::string DefaultGravityName()

Definition GravityModel.cpp:369

Math::real GeoidHeight(real lat, real lon) const

Definition GravityModel.cpp:290

Math::real T(real X, real Y, real Z, real &deltaX, real &deltaY, real &deltaZ) const

Math::real Gravity(real lat, real lon, real h, real &gx, real &gy, real &gz) const

Definition GravityModel.cpp:304

Math::real V(real X, real Y, real Z, real &GX, real &GY, real &GZ) const

Definition GravityModel.cpp:243

Math::real W(real X, real Y, real Z, real &gX, real &gY, real &gZ) const

Definition GravityModel.cpp:255

Mathematical functions needed by GeographicLib.

const Geocentric & Earth() const

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

Math::real MassConstant() const

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

Math::real SurfaceGravity(real lat) const

static void readcoeffs(std::istream &stream, int &N, int &M, std::vector< real > &C, std::vector< real > &S, bool truncate=false)

Spherical harmonic series with a correction to the coefficients.

CircularEngine Circle(real tau, real p, real z, bool gradp) const

Spherical harmonic series.

const SphericalEngine::coeff & Coefficients() const

CircularEngine Circle(real p, real z, bool gradp) const

static bool ParseLine(const std::string &line, std::string &key, std::string &value, char equals='\0', char comment='#')

Namespace for GeographicLib.