00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <GeographicLib/GravityCircle.hpp>
00011 #include <fstream>
00012 #include <sstream>
00013 #include <GeographicLib/Geocentric.hpp>
00014
00015 namespace GeographicLib {
00016
00017 using namespace std;
00018
00019 Math::real GravityCircle::Gravity(real lon, real& gx, real& gy, real& gz)
00020 const {
00021 real clam, slam, M[Geocentric::dim2_];
00022 CircularEngine::cossin(lon, clam, slam);
00023 real Wres = W(clam, slam, gx, gy, gz);
00024 Geocentric::Rotation(_sphi, _cphi, slam, clam, M);
00025 Geocentric::Unrotate(M, gx, gy, gz, gx, gy, gz);
00026 return Wres;
00027 }
00028
00029 Math::real GravityCircle::Disturbance(real lon, real& deltax, real& deltay,
00030 real& deltaz) const {
00031 real clam, slam, M[Geocentric::dim2_];
00032 CircularEngine::cossin(lon, clam, slam);
00033 real Tres = InternalT(clam, slam, deltax, deltay, deltaz, true, true);
00034 Geocentric::Rotation(_sphi, _cphi, slam, clam, M);
00035 Geocentric::Unrotate(M, deltax, deltay, deltaz, deltax, deltay, deltaz);
00036 return Tres;
00037 }
00038
00039 Math::real GravityCircle::GeoidHeight(real lon) const {
00040 if ((_caps & GEOID_HEIGHT) != GEOID_HEIGHT)
00041 return Math::NaN();
00042 real clam, slam, dummy;
00043 CircularEngine::cossin(lon, clam, slam);
00044 real T = InternalT(clam, slam, dummy, dummy, dummy, false, false);
00045 real correction = _corrmult * _correction(clam, slam);
00046 return T/_gamma0 + correction;
00047 }
00048
00049 void GravityCircle::SphericalAnomaly(real lon,
00050 real& Dg01, real& xi, real& eta)
00051 const {
00052 if ((_caps & SPHERICAL_ANOMALY) != SPHERICAL_ANOMALY) {
00053 Dg01 = xi = eta = Math::NaN();
00054 return;
00055 }
00056 real clam, slam;
00057 CircularEngine::cossin(lon, clam, slam);
00058 real
00059 deltax, deltay, deltaz,
00060 T = InternalT(clam, slam, deltax, deltay, deltaz, true, false);
00061
00062 real MC[Geocentric::dim2_];
00063 Geocentric::Rotation(_spsi, _cpsi, slam, clam, MC);
00064 Geocentric::Unrotate(MC, deltax, deltay, deltaz, deltax, deltay, deltaz);
00065
00066 Dg01 = - deltaz - 2 * T * _invR;
00067 xi = -(deltay/_gamma) / Math::degree();
00068 eta = -(deltax/_gamma) / Math::degree();
00069 }
00070
00071 Math::real GravityCircle::W(real clam, real slam,
00072 real& gX, real& gY, real& gZ) const {
00073 real Wres = V(clam, slam, gX, gY, gZ) + _frot * _Px / 2;
00074 gX += _frot * clam;
00075 gY += _frot * slam;
00076 return Wres;
00077 }
00078
00079 Math::real GravityCircle::V(real clam, real slam,
00080 real& GX, real& GY, real& GZ)
00081 const {
00082 if ((_caps & GRAVITY) != GRAVITY) {
00083 GX = GY = GZ = Math::NaN();
00084 return Math::NaN();
00085 }
00086 real
00087 Vres = _gravitational(clam, slam, GX, GY, GZ),
00088 f = _GMmodel / _amodel;
00089 Vres *= f;
00090 GX *= f;
00091 GY *= f;
00092 GZ *= f;
00093 return Vres;
00094 }
00095
00096 Math::real GravityCircle::InternalT(real clam, real slam,
00097 real& deltaX, real& deltaY, real& deltaZ,
00098 bool gradp, bool correct) const {
00099 if (gradp) {
00100 if ((_caps & DISTURBANCE) != DISTURBANCE) {
00101 deltaX = deltaY = deltaZ = Math::NaN();
00102 return Math::NaN();
00103 }
00104 } else {
00105 if ((_caps & DISTURBING_POTENTIAL) != DISTURBING_POTENTIAL)
00106 return Math::NaN();
00107 }
00108 if (_dzonal0 == 0)
00109 correct = false;
00110 real T = (gradp
00111 ? _disturbing(clam, slam, deltaX, deltaY, deltaZ)
00112 : _disturbing(clam, slam));
00113 T = (T / _amodel - (correct ? _dzonal0 : 0) * _invR) * _GMmodel;
00114 if (gradp) {
00115 real f = _GMmodel / _amodel;
00116 deltaX *= f;
00117 deltaY *= f;
00118 deltaZ *= f;
00119 if (correct) {
00120 real r3 = _GMmodel * _dzonal0 * _invR * _invR * _invR;
00121 deltaX += _Px * clam * r3;
00122 deltaY += _Px * slam * r3;
00123 deltaZ += _Z * r3;
00124 }
00125 }
00126 return T;
00127 }
00128
00129 }