GeographicLib  1.51
GravityCircle.cpp
Go to the documentation of this file.
1 /**
2  * \file GravityCircle.cpp
3  * \brief Implementation for GeographicLib::GravityCircle class
4  *
5  * Copyright (c) Charles Karney (2011-2020) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
11 #include <fstream>
12 #include <sstream>
14 
15 namespace GeographicLib {
16 
17  using namespace std;
18 
19  GravityCircle::GravityCircle(mask caps, real a, real f, real lat, real h,
20  real Z, real P, real cphi, real sphi,
21  real amodel, real GMmodel,
22  real dzonal0, real corrmult,
23  real gamma0, real gamma, real frot,
24  const CircularEngine& gravitational,
25  const CircularEngine& disturbing,
26  const CircularEngine& correction)
27  : _caps(caps)
28  , _a(a)
29  , _f(f)
30  , _lat(Math::LatFix(lat))
31  , _h(h)
32  , _Z(Z)
33  , _Px(P)
34  , _invR(1 / hypot(_Px, _Z))
35  , _cpsi(_Px * _invR)
36  , _spsi(_Z * _invR)
37  , _cphi(cphi)
38  , _sphi(sphi)
39  , _amodel(amodel)
40  , _GMmodel(GMmodel)
41  , _dzonal0(dzonal0)
42  , _corrmult(corrmult)
43  , _gamma0(gamma0)
44  , _gamma(gamma)
45  , _frot(frot)
46  , _gravitational(gravitational)
47  , _disturbing(disturbing)
48  , _correction(correction)
49  {}
50 
51  Math::real GravityCircle::Gravity(real lon,
52  real& gx, real& gy, real& gz) const {
53  real slam, clam, M[Geocentric::dim2_];
54  Math::sincosd(lon, slam, clam);
55  real Wres = W(slam, clam, gx, gy, gz);
56  Geocentric::Rotation(_sphi, _cphi, slam, clam, M);
57  Geocentric::Unrotate(M, gx, gy, gz, gx, gy, gz);
58  return Wres;
59  }
60 
61  Math::real GravityCircle::Disturbance(real lon, real& deltax, real& deltay,
62  real& deltaz) const {
63  real slam, clam, M[Geocentric::dim2_];
64  Math::sincosd(lon, slam, clam);
65  real Tres = InternalT(slam, clam, deltax, deltay, deltaz, true, true);
66  Geocentric::Rotation(_sphi, _cphi, slam, clam, M);
67  Geocentric::Unrotate(M, deltax, deltay, deltaz, deltax, deltay, deltaz);
68  return Tres;
69  }
70 
71  Math::real GravityCircle::GeoidHeight(real lon) const {
72  if ((_caps & GEOID_HEIGHT) != GEOID_HEIGHT)
73  return Math::NaN();
74  real slam, clam, dummy;
75  Math::sincosd(lon, slam, clam);
76  real T = InternalT(slam, clam, dummy, dummy, dummy, false, false);
77  real correction = _corrmult * _correction(slam, clam);
78  return T/_gamma0 + correction;
79  }
80 
81  void GravityCircle::SphericalAnomaly(real lon,
82  real& Dg01, real& xi, real& eta) const {
83  if ((_caps & SPHERICAL_ANOMALY) != SPHERICAL_ANOMALY) {
84  Dg01 = xi = eta = Math::NaN();
85  return;
86  }
87  real slam, clam;
88  Math::sincosd(lon, slam, clam);
89  real
90  deltax, deltay, deltaz,
91  T = InternalT(slam, clam, deltax, deltay, deltaz, true, false);
92  // Rotate cartesian into spherical coordinates
93  real MC[Geocentric::dim2_];
94  Geocentric::Rotation(_spsi, _cpsi, slam, clam, MC);
95  Geocentric::Unrotate(MC, deltax, deltay, deltaz, deltax, deltay, deltaz);
96  // H+M, Eq 2-151c
97  Dg01 = - deltaz - 2 * T * _invR;
98  xi = -(deltay/_gamma) / Math::degree();
99  eta = -(deltax/_gamma) / Math::degree();
100  }
101 
102  Math::real GravityCircle::W(real slam, real clam,
103  real& gX, real& gY, real& gZ) const {
104  real Wres = V(slam, clam, gX, gY, gZ) + _frot * _Px / 2;
105  gX += _frot * clam;
106  gY += _frot * slam;
107  return Wres;
108  }
109 
110  Math::real GravityCircle::V(real slam, real clam,
111  real& GX, real& GY, real& GZ) const {
112  if ((_caps & GRAVITY) != GRAVITY) {
113  GX = GY = GZ = Math::NaN();
114  return Math::NaN();
115  }
116  real
117  Vres = _gravitational(slam, clam, GX, GY, GZ),
118  f = _GMmodel / _amodel;
119  Vres *= f;
120  GX *= f;
121  GY *= f;
122  GZ *= f;
123  return Vres;
124  }
125 
126  Math::real GravityCircle::InternalT(real slam, real clam,
127  real& deltaX, real& deltaY, real& deltaZ,
128  bool gradp, bool correct) const {
129  if (gradp) {
130  if ((_caps & DISTURBANCE) != DISTURBANCE) {
131  deltaX = deltaY = deltaZ = Math::NaN();
132  return Math::NaN();
133  }
134  } else {
135  if ((_caps & DISTURBING_POTENTIAL) != DISTURBING_POTENTIAL)
136  return Math::NaN();
137  }
138  if (_dzonal0 == 0)
139  correct = false;
140  real T = (gradp
141  ? _disturbing(slam, clam, deltaX, deltaY, deltaZ)
142  : _disturbing(slam, clam));
143  T = (T / _amodel - (correct ? _dzonal0 : 0) * _invR) * _GMmodel;
144  if (gradp) {
145  real f = _GMmodel / _amodel;
146  deltaX *= f;
147  deltaY *= f;
148  deltaZ *= f;
149  if (correct) {
150  real r3 = _GMmodel * _dzonal0 * _invR * _invR * _invR;
151  deltaX += _Px * clam * r3;
152  deltaY += _Px * slam * r3;
153  deltaZ += _Z * r3;
154  }
155  }
156  return T;
157  }
158 
159 } // namespace GeographicLib
real
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
GeographicLib
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
GeographicLib::Math::real
double real
Definition: Math.hpp:99
GravityCircle.hpp
Header for GeographicLib::GravityCircle class.
Geocentric.hpp
Header for GeographicLib::Geocentric class.
std
Definition: NearestNeighbor.hpp:813
GeographicLib::GravityCircle::GravityCircle
GravityCircle()
Definition: GravityCircle.hpp:83