GeographicLib  1.21
GravityCircle.cpp
Go to the documentation of this file.
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