src/Environment/gravity.cxx

   1 // gravity.cxx -- interface for earth gravitational model
   2 //
   3 // Written by Torsten Dreyer, June 2011
   4 //
   5 // Copyright (C) 2011  Torsten Dreyer - torsten (at) t3r _dot_ de
   6 //
   7 // This program is free software; you can redistribute it and/or
   8 // modify it under the terms of the GNU General Public License as
   9 // published by the Free Software Foundation; either version 2 of the
  10 // License, or (at your option) any later version.
  11 //
  12 // This program is distributed in the hope that it will be useful, but
  13 // WITHOUT ANY WARRANTY; without even the implied warranty of
  14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15 // General Public License for more details.
  16 //
  17 // You should have received a copy of the GNU General Public License
  18 // along with this program; if not, write to the Free Software
  19 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  20 //
  21
  22 #include "gravity.hxx"
  23
  24 #include <simgear/structure/exception.hxx>
  25
  26 namespace Environment {
  27
  28 /*
  29 http://de.wikipedia.org/wiki/Normalschwereformel
  30 */
  31 class Somigliana : public Gravity {
  32 public:
  33     Somigliana();
  34     virtual ~Somigliana();
  35     virtual double getGravity( const SGGeod & position ) const;
  36 };
  37
  38 Somigliana::Somigliana()
  39 {
  40 }
  41
  42 Somigliana::~Somigliana()
  43 {
  44 }
  45
  46 double Somigliana::getGravity( const SGGeod & position ) const
  47 {
  48 // Geodetic Reference System 1980 parameter
  49 #define A 6378137.0 // equatorial radius of earth
  50 #define B 6356752.3141 // semiminor axis
  51 #define AGA (A*9.7803267715) // A times normal gravity at equator
  52 #define BGB (B*9.8321863685) // B times normal gravity at pole
  53   // forumla of Somigliana
  54   double cosphi = ::cos(position.getLatitudeRad());
  55   double cos2phi = cosphi*cosphi;
  56   double sinphi = ::sin(position.getLatitudeRad());
  57   double sin2phi = sinphi*sinphi;
  58   double g0 = (AGA * cos2phi + BGB * sin2phi) / sqrt( A*A*cos2phi+B*B*sin2phi );
  59
  60   static const double k1 = 3.15704e-7;
  61   static const double k2 = 2.10269e-9;
  62   static const double k3 = 7.37452e-14;
  63
  64   double h = position.getElevationM();
  65
  66   return g0*(1-(k1-k2*sin2phi)*h+k3*h*h);
  67 }
  68
  69 static Somigliana _somigliana;
  70
  71 /* --------------------- Gravity implementation --------------------- */
  72 Gravity * Gravity::_instance = NULL;
  73
  74 Gravity::~Gravity()
  75 {
  76 }
  77
  78 //double Gravity::getGravity( const SGGeoc & position ) = 0;
  79
  80 const Gravity * Gravity::instance()
  81 {
  82     if( _instance == NULL )
  83         _instance = &_somigliana;
  84
  85     return _instance;
  86 }
  87
  88 } // namespace