Compute gravitational acceleration

Compute gravitational acceleration based on the formula of
Somigliana corrected for altitude and provide the value
in property /environment/gravitational-acceleration-mps2

projects/VC90/FlightGear/FlightGear.vcproj

@@ -3437,6 +3437,14 @@
 				RelativePath="..\..\..\src\Environment\presets.hxx"
 				>
 			</File>
+			<File
+				RelativePath="..\..\..\src\Environment\gravity.cxx"
+				>
+			</File>
+			<File
+				RelativePath="..\..\..\src\Environment\gravity.hxx"
+				>
+			</File>
 		</Filter>
 		<Filter
 			Name="Lib_Model"

src/Environment/CMakeLists.txt

@@ -16,6 +16,7 @@ set(SOURCES
 	ridge_lift.cxx
 	terrainsampler.cxx
 	presets.cxx
+	gravity.cxx
 	)
 	
 flightgear_component(Environment "${SOURCES}")

src/Environment/Makefile.am

@@ -17,6 +17,7 @@ libEnvironment_a_SOURCES = \
 	ridge_lift.cxx ridge_lift.hxx \
 	ephemeris.cxx ephemeris.hxx \
 	terrainsampler.cxx terrainsampler.cxx \
-	presets.cxx presets.hxx
+	presets.cxx presets.hxx \
+	gravity.cxx gravity.hxx
 
 INCLUDES = -I$(top_srcdir) -I$(top_srcdir)/src

src/Environment/environment_mgr.cxx

@@ -42,6 +42,7 @@
 #include "ridge_lift.hxx"
 #include "terrainsampler.hxx"
 #include "Airports/simple.hxx"
+#include "gravity.hxx"
 
 class SGSky;
 extern SGSky *thesky;
@@ -50,7 +51,7 @@ FGEnvironmentMgr::FGEnvironmentMgr () :
   _environment(new FGEnvironment()),
   fgClouds(new FGClouds()),
   _cloudLayersDirty(true),
-  _altitudeNode(fgGetNode("/position/altitude-ft", true)),
+  _altitude_n(fgGetNode("/position/altitude-ft", true)),
   _longitude_n(fgGetNode( "/position/longitude-deg", true )),
   _latitude_n( fgGetNode( "/position/latitude-deg", true )),
   _positionTimeToLive(0.0)
@@ -98,10 +99,11 @@ FGEnvironmentMgr::init ()
   SGSubsystemGroup::init();
   fgClouds->Init();
 
+  // FIXME: is this really part of the environment_mgr?
   // Initialize the longitude, latitude and altitude to the initial position
   // of the aircraft so that the atmospheric properties (pressure, temperature
   // and density) can be initialized accordingly.
-  _altitudeNode->setDoubleValue(fgGetDouble("/sim/presets/altitude-ft"));
+  _altitude_n->setDoubleValue(fgGetDouble("/sim/presets/altitude-ft"));
   _longitude_n->setDoubleValue(fgGetDouble("/sim/presets/longitude-deg"));
   _latitude_n->setDoubleValue(fgGetDouble("/sim/presets/latitude-deg"));
 }
@@ -201,7 +203,7 @@ FGEnvironmentMgr::update (double dt)
 {
   SGSubsystemGroup::update(dt);
 
-  _environment->set_elevation_ft( _altitudeNode->getDoubleValue() );
+  _environment->set_elevation_ft( _altitude_n->getDoubleValue() );
 
   simgear::Particles::setWindFrom( _environment->get_wind_from_heading_deg(),
                                    _environment->get_wind_speed_kt() );
@@ -210,6 +212,14 @@ FGEnvironmentMgr::update (double dt)
     fgClouds->set_update_event( fgClouds->get_update_event()+1 );
   }
 
+
+  fgSetDouble( "/environment/gravitational-acceleration-mps2", 
+    Environment::Gravity::instance()->getGravity(SGGeod::fromDegFt(
+      _longitude_n->getDoubleValue(),
+      _latitude_n->getDoubleValue(),
+      _altitude_n->getDoubleValue()
+  )));
+
   _positionTimeToLive -= dt;
   if( _positionTimeToLive <= 0.0 )
   {

src/Environment/environment_mgr.hxx

@@ -95,7 +95,7 @@ private:
   FGEnvironment * _environment;	// always the same, for now
   FGClouds *fgClouds;
   bool _cloudLayersDirty;
-  SGPropertyNode_ptr _altitudeNode;
+  SGPropertyNode_ptr _altitude_n;
   SGPropertyNode_ptr _longitude_n;
   SGPropertyNode_ptr _latitude_n;
   double _positionTimeToLive;

src/Environment/gravity.cxx

@@ -0,0 +1,89 @@
+// gravity.cxx -- interface for earth gravitational model
+//
+// Written by Torsten Dreyer, June 2011
+//
+// Copyright (C) 2011  Torsten Dreyer - torsten (at) t3r _dot_ de
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of the
+// License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+
+#include "gravity.hxx"
+
+#include <simgear/structure/exception.hxx>
+
+namespace Environment {
+
+/*
+http://de.wikipedia.org/wiki/Normalschwereformel
+*/
+class Somigliana : public Gravity {
+public:
+    Somigliana();
+    virtual ~Somigliana();
+    virtual double getGravity( const SGGeod & position ) const;
+};
+
+Somigliana::Somigliana()
+{
+}
+
+Somigliana::~Somigliana()
+{
+}
+#include <stdio.h>
+
+double Somigliana::getGravity( const SGGeod & position ) const
+{
+// Geodetic Reference System 1980 parameter
+#define A 6378137.0 // equatorial radius of earth
+#define B 6356752.3141 // semiminor axis
+#define AGA (A*9.7803267715) // A times normal gravity at equator
+#define BGB (B*9.8321863685) // B times normal gravity at pole
+  // forumla of Somigliana
+  double cosphi = ::cos(position.getLatitudeRad());
+  double cos2phi = cosphi*cosphi;
+  double sinphi = ::sin(position.getLatitudeRad());
+  double sin2phi = sinphi*sinphi;
+  double g0 = (AGA * cos2phi + BGB * sin2phi) / sqrt( A*A*cos2phi+B*B*sin2phi );
+
+  static const double k1 = 3.15704e-7;
+  static const double k2 = 2.10269e-9;
+  static const double k3 = 7.37452e-14;
+
+  double h = position.getElevationM();
+  
+  return g0*(1-(k1-k2*sin2phi)*h+k3*h*h);
+}
+
+static Somigliana _somigliana;
+
+/* --------------------- Gravity implementation --------------------- */
+Gravity * Gravity::_instance = NULL;
+
+Gravity::~Gravity()
+{
+}
+
+//double Gravity::getGravity( const SGGeoc & position ) = 0;
+
+const Gravity * Gravity::instance()
+{
+    if( _instance == NULL )
+        _instance = &_somigliana;
+
+    return _instance;
+}
+    
+} // namespace

src/Environment/gravity.hxx

@@ -0,0 +1,43 @@
+// gravity.hxx -- interface for earth gravitational model
+//
+// Written by Torsten Dreyer, June 2011
+//
+// Copyright (C) 2011  Torsten Dreyer - torsten (at) t3r _dot_ de
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of the
+// License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+
+#ifndef __GRAVITY_HXX
+#define __GRAVITY_HXX
+
+#include <simgear/math/SGMath.hxx>
+
+namespace Environment {
+
+class Gravity 
+{
+public:
+    virtual ~Gravity();
+    virtual double getGravity( const SGGeod & position ) const = 0;
+
+    const static Gravity * instance();
+
+private:
+    static Gravity * _instance;
+    
+};
+
+} // namespace
+#endif // __GRAVITY_HXX