GeographicLib
1.21
|
00001 /** 00002 * \file GravityCircle.cpp 00003 * \brief Implementation for GeographicLib::GravityCircle class 00004 * 00005 * Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed under 00006 * the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #include <GeographicLib/GravityCircle.hpp> 00011 #include <fstream> 00012 #include <sstream> 00013 #include <GeographicLib/Geocentric.hpp> 00014 00015 #define GEOGRAPHICLIB_GRAVITYCIRCLE_CPP \ 00016 "$Id: 94c2bce8fa9f379b5e3d0e176c89b4989762ef3a $" 00017 00018 RCSID_DECL(GEOGRAPHICLIB_GRAVITYCIRCLE_CPP) 00019 RCSID_DECL(GEOGRAPHICLIB_GRAVITYCIRCLE_HPP) 00020 00021 #define GRAVITY_DEFAULT_PATH "/home/ckarney/geographiclib/gravity" 00022 00023 namespace GeographicLib { 00024 00025 using namespace std; 00026 00027 Math::real GravityCircle::Gravity(real lon, real& gx, real& gy, real& gz) 00028 const throw() { 00029 real clam, slam, M[Geocentric::dim2_]; 00030 CircularEngine::cossin(lon, clam, slam); 00031 real Wres = W(clam, slam, gx, gy, gz); 00032 Geocentric::Rotation(_sphi, _cphi, slam, clam, M); 00033 Geocentric::Unrotate(M, gx, gy, gz, gx, gy, gz); 00034 return Wres; 00035 } 00036 00037 Math::real GravityCircle::Disturbance(real lon, real& deltax, real& deltay, 00038 real& deltaz) const throw() { 00039 real clam, slam, M[Geocentric::dim2_]; 00040 CircularEngine::cossin(lon, clam, slam); 00041 real Tres = InternalT(clam, slam, deltax, deltay, deltaz, true, true); 00042 Geocentric::Rotation(_sphi, _cphi, slam, clam, M); 00043 Geocentric::Unrotate(M, deltax, deltay, deltaz, deltax, deltay, deltaz); 00044 return Tres; 00045 } 00046 00047 Math::real GravityCircle::GeoidHeight(real lon) const throw() { 00048 if ((_caps & GEOID_HEIGHT) != GEOID_HEIGHT) 00049 return Math::NaN<real>(); 00050 real clam, slam, dummy; 00051 CircularEngine::cossin(lon, clam, slam); 00052 real T = InternalT(clam, slam, dummy, dummy, dummy, false, false); 00053 real correction = _corrmult * _correction(clam, slam); 00054 return T/_gamma0 + correction; 00055 } 00056 00057 void GravityCircle::SphericalAnomaly(real lon, 00058 real& Dg01, real& xi, real& eta) 00059 const throw() { 00060 if ((_caps & SPHERICAL_ANOMALY) != SPHERICAL_ANOMALY) { 00061 Dg01 = xi = eta = Math::NaN<real>(); 00062 return; 00063 } 00064 real clam, slam; 00065 CircularEngine::cossin(lon, clam, slam); 00066 real 00067 deltax, deltay, deltaz, 00068 T = InternalT(clam, slam, deltax, deltay, deltaz, true, false); 00069 // Rotate cartesian into spherical coordinates 00070 real MC[Geocentric::dim2_]; 00071 Geocentric::Rotation(_spsi, _cpsi, slam, clam, MC); 00072 Geocentric::Unrotate(MC, deltax, deltay, deltaz, deltax, deltay, deltaz); 00073 // H+M, Eq 2-151c 00074 Dg01 = - deltaz - 2 * T * _invR; 00075 xi = -(deltay/_gamma) / Math::degree<real>(); 00076 eta = -(deltax/_gamma) / Math::degree<real>(); 00077 } 00078 00079 Math::real GravityCircle::W(real clam, real slam, 00080 real& gX, real& gY, real& gZ) const throw() { 00081 real Wres = V(clam, slam, gX, gY, gZ) + _frot * _P / 2; 00082 gX += _frot * clam; 00083 gY += _frot * slam; 00084 return Wres; 00085 } 00086 00087 Math::real GravityCircle::V(real clam, real slam, 00088 real& GX, real& GY, real& GZ) 00089 const throw() { 00090 if ((_caps & GRAVITY) != GRAVITY) { 00091 GX = GY = GZ = Math::NaN<real>(); 00092 return Math::NaN<real>(); 00093 } 00094 real 00095 Vres = _gravitational(clam, slam, GX, GY, GZ), 00096 f = _GMmodel / _amodel; 00097 Vres *= f; 00098 GX *= f; 00099 GY *= f; 00100 GZ *= f; 00101 return Vres; 00102 } 00103 00104 Math::real GravityCircle::InternalT(real clam, real slam, 00105 real& deltaX, real& deltaY, real& deltaZ, 00106 bool gradp, bool correct) const throw() { 00107 if (gradp) { 00108 if ((_caps & DISTURBANCE) != DISTURBANCE) { 00109 deltaX = deltaY = deltaZ = Math::NaN<real>(); 00110 return Math::NaN<real>(); 00111 } 00112 } else { 00113 if ((_caps & DISTURBING_POTENTIAL) != DISTURBING_POTENTIAL) 00114 return Math::NaN<real>(); 00115 } 00116 if (_dzonal0 == 0) 00117 correct = false; 00118 real T = (gradp 00119 ? _disturbing(clam, slam, deltaX, deltaY, deltaZ) 00120 : _disturbing(clam, slam)); 00121 T = (T / _amodel - (correct ? _dzonal0 : 0) * _invR) * _GMmodel; 00122 if (gradp) { 00123 real f = _GMmodel / _amodel; 00124 deltaX *= f; 00125 deltaY *= f; 00126 deltaZ *= f; 00127 if (correct) { 00128 real r3 = _GMmodel * _dzonal0 * _invR * _invR * _invR; 00129 deltaX += _P * clam * r3; 00130 deltaY += _P * slam * r3; 00131 deltaZ += _Z * r3; 00132 } 00133 } 00134 return T; 00135 } 00136 00137 } // namespace GeographicLib