1 // gravity.cxx -- interface for earth gravitational model
3 // Written by Torsten Dreyer, June 2011
5 // Copyright (C) 2011 Torsten Dreyer - torsten (at) t3r _dot_ de
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.
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.
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.
22 #include "gravity.hxx"
24 #include <simgear/structure/exception.hxx>
26 namespace Environment {
29 http://de.wikipedia.org/wiki/Normalschwereformel
31 class Somigliana : public Gravity {
34 virtual ~Somigliana();
35 virtual double getGravity( const SGGeod & position ) const;
38 Somigliana::Somigliana()
42 Somigliana::~Somigliana()
46 double Somigliana::getGravity( const SGGeod & position ) const
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 );
60 static const double k1 = 3.15704e-7;
61 static const double k2 = 2.10269e-9;
62 static const double k3 = 7.37452e-14;
64 double h = position.getElevationM();
66 return g0*(1-(k1-k2*sin2phi)*h+k3*h*h);
69 static Somigliana _somigliana;
71 /* --------------------- Gravity implementation --------------------- */
72 Gravity * Gravity::_instance = NULL;
78 //double Gravity::getGravity( const SGGeoc & position ) = 0;
80 const Gravity * Gravity::instance()
82 if( _instance == NULL )
83 _instance = &_somigliana;