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 #include <stdio.h>
  46
  47 double Somigliana::getGravity( const SGGeod & position ) const
  48 {
  49 // Geodetic Reference System 1980 parameter
  50 #define A 6378137.0 // equatorial radius of earth
  51 #define B 6356752.3141 // semiminor axis
  52 #define AGA (A*9.7803267715) // A times normal gravity at equator
  53 #define BGB (B*9.8321863685) // B times normal gravity at pole
  54   // forumla of Somigliana
  55   double cosphi = ::cos(position.getLatitudeRad());
  56   double cos2phi = cosphi*cosphi;
  57   double sinphi = ::sin(position.getLatitudeRad());
  58   double sin2phi = sinphi*sinphi;
  59   double g0 = (AGA * cos2phi + BGB * sin2phi) / sqrt( A*A*cos2phi+B*B*sin2phi );
  60
  61   static const double k1 = 3.15704e-7;
  62   static const double k2 = 2.10269e-9;
  63   static const double k3 = 7.37452e-14;
  64
  65   double h = position.getElevationM();
  66
  67   return g0*(1-(k1-k2*sin2phi)*h+k3*h*h);
  68 }
  69
  70 static Somigliana _somigliana;
  71
  72 /* --------------------- Gravity implementation --------------------- */
  73 Gravity * Gravity::_instance = NULL;
  74
  75 Gravity::~Gravity()
  76 {
  77 }
  78
  79 //double Gravity::getGravity( const SGGeoc & position ) = 0;
  80
  81 const Gravity * Gravity::instance()
  82 {
  83     if( _instance == NULL )
  84         _instance = &_somigliana;
  85
  86     return _instance;
  87 }
  88
  89 } // namespace