From: Eric van den Berg <evan_den_berg@hotmail.com>
Date: Mon, 26 Nov 2012 17:59:48 +0000 (+0100)
Subject: calculate /systems/pitot/total-pressure on the bases of /velocities/mach
X-Git-Url: https://git.mxchange.org/?a=commitdiff_plain;h=5a8e6aeb8626db27fc606da1e0c8dfa913c88265;p=flightgear.git

calculate /systems/pitot/total-pressure on the bases of /velocities/mach
added measured-total-pressure property to account for a normal
shockwave in front of pitot tube at supersonic speeds. momentarily not
used
small unit conversion correction of inHg to hPa in altimeter
Base airspeed calculations on impact pressure in instrumentation
/airspeed_indicator
---

diff --git a/src/Instrumentation/airspeed_indicator.cxx b/src/Instrumentation/airspeed_indicator.cxx
index 223c4d774..3c339451c 100644
--- a/src/Instrumentation/airspeed_indicator.cxx
+++ b/src/Instrumentation/airspeed_indicator.cxx
@@ -1,5 +1,6 @@
 // airspeed_indicator.cxx - a regular pitot-static airspeed indicator.
 // Written by David Megginson, started 2002.
+// Last modified by Eric van den Berg, 09 Dec 2012
 //
 // This file is in the Public Domain and comes with no warranty.
 
@@ -9,11 +10,14 @@
 
 #include <math.h>
 
+#include <simgear/constants.h>
 #include <simgear/math/interpolater.hxx>
 
 #include "airspeed_indicator.hxx"
 #include <Main/fg_props.hxx>
 #include <Main/util.hxx>
+#include <Environment/environment_mgr.hxx>
+#include <Environment/environment.hxx>
 
 
 // A higher number means more responsive.
@@ -31,6 +35,7 @@ AirspeedIndicator::AirspeedIndicator ( SGPropertyNode *node )
     _mach_limit(node->getDoubleValue("mach-limit", 0.48)),
     _alt_threshold(node->getDoubleValue("alt-threshold", 13200))
 {
+  _environmentManager = NULL;
 }
 
 AirspeedIndicator::~AirspeedIndicator ()
@@ -48,10 +53,6 @@ AirspeedIndicator::init ()
     _total_pressure_node = fgGetNode(_total_pressure.c_str(), true);
     _static_pressure_node = fgGetNode(_static_pressure.c_str(), true);
     _density_node = fgGetNode("/environment/density-slugft3", true);
-    
-    _sea_level_pressure_node = fgGetNode("/environment/pressure-sea-level-inhg", true);
-    _oat_celsius_node = fgGetNode("/environment/temperature-degc", true);
-  
     _speed_node = node->getChild("indicated-speed-kt", 0, true);
     _tas_node = node->getChild("true-speed-kt", 0, true);
     _mach_node = node->getChild("indicated-mach", 0, true);
@@ -77,6 +78,8 @@ AirspeedIndicator::init ()
         _airspeed_limit = node->getChild("airspeed-limit-kt", 0, true);
         _pressure_alt = fgGetNode(_pressure_alt_source.c_str(), true);
     }
+    
+    _environmentManager = (FGEnvironmentMgr*) globals->get_subsystem("environment");
 }
 
 void
@@ -85,14 +88,6 @@ AirspeedIndicator::reinit ()
     _speed_node->setDoubleValue(0.0);
 }
 
-#ifndef FPSTOKTS
-# define FPSTOKTS 0.592484
-#endif
-
-#ifndef INHGTOPSF
-# define INHGTOPSF (2116.217/29.9212)
-#endif
-
 void
 AirspeedIndicator::update (double dt)
 {
@@ -100,24 +95,25 @@ AirspeedIndicator::update (double dt)
         return;
     }
     
-    double pt = _total_pressure_node->getDoubleValue() * INHGTOPSF;
-    double p = _static_pressure_node->getDoubleValue() * INHGTOPSF;
-    double r = _density_node->getDoubleValue();
-    double q = ( pt - p );  // dynamic pressure
+    double pt = _total_pressure_node->getDoubleValue() ;
+    double p = _static_pressure_node->getDoubleValue() ;
+    double qc = ( pt - p ) * SG_INHG_TO_PA ;  // Impact pressure in Pa, _not_ to be confused with dynamic pressure!!!
 
-    // Now, reverse the equation (normalize dynamic pressure to
+    // Now, reverse the equation (normalize impact pressure to
     // avoid "nan" results from sqrt)
-    if ( q < 0 ) { q = 0.0; }
-    double v_fps = sqrt((2 * q) / r);
+    qc = std::max( qc , 0.0 );
+    // Calibrated airspeed (using compressible aerodynamics) based on impact pressure qc in m/s
+    // Using calibrated airspeed as indicated airspeed, neglecting any airspeed indicator errors.
+    double v_cal = sqrt( 7 * SG_p0_Pa/SG_rho0_kg_p_m3 * ( pow( 1 + qc/SG_p0_Pa  , 1/3.5 )  -1 ) );
 
                             // Publish the indicated airspeed
     double last_speed_kt = _speed_node->getDoubleValue();
-    double current_speed_kt = v_fps * FPSTOKTS;
+    double current_speed_kt = v_cal * SG_MPS_TO_KT;
     double filtered_speed = fgGetLowPass(last_speed_kt,
                                              current_speed_kt,
                                              dt * RESPONSIVENESS);
     _speed_node->setDoubleValue(filtered_speed);
-    computeMach(filtered_speed);
+    computeMach();
 
     if (!_has_overspeed) {
         return;
@@ -133,24 +129,32 @@ AirspeedIndicator::update (double dt)
 }
 
 void
-AirspeedIndicator::computeMach(double ias)
+AirspeedIndicator::computeMach()
 {
-   
-  // derived from http://williams.best.vwh.net/avform.htm#Mach
-  // names here are picked to be consistent with those formulae!
-
-  double oatK = _oat_celsius_node->getDoubleValue() + 273.15; // OAT in Kelvin
-  double CS = 38.967854 * sqrt(oatK); // speed-of-sound in knots at altitude
-  double CS_0 = 661.4786; // speed-of-sound in knots at sea-level
-  double P_0 = _sea_level_pressure_node->getDoubleValue();
-  double P = _static_pressure_node->getDoubleValue();
+  if (!_environmentManager) {
+    return;
+  }
+  
+  FGEnvironment env(_environmentManager->getEnvironment());
+  
+  double oatK = env.get_temperature_degc() + SG_T0_K - 15.0 ;         // OAT in Kelvin
+  oatK = std::max( oatK , 0.001 );                                // should never happen, but just in case someone flies into space...
+  double c = sqrt(SG_gamma * SG_R_m2_p_s2_p_K * oatK);                 // speed-of-sound in m/s at aircraft position
+  double pt = _total_pressure_node->getDoubleValue() * SG_INHG_TO_PA;  // total pressure in Pa
+  double p = _static_pressure_node->getDoubleValue() * SG_INHG_TO_PA;  // static pressure in Pa
+  p = std::max( p , 0.001 );                                           // should never happen, but just in case someone flies into space...
+  double rho = _density_node->getDoubleValue() * SG_SLUGFT3_TO_KGPM3;  // air density in kg/m3
+  rho = std::max( rho , 0.001 );                                      // should never happen, but just in case someone flies into space...
   
-  double DP = P_0 * (pow(1 + 0.2*pow(ias/CS_0, 2), 3.5) - 1);
-  double mach = pow(5 * ( pow(DP/P + 1, 2.0/7.0) -1) , 0.5);
+  // true airspeed in m/s
+  pt = std::max( pt , p );
+  double V_true = sqrt( 7 * p/rho * (pow( 1 + (pt-p)/p , 1/3.5 ) -1 ) );
+  // Mach number; _see notes in systems/pitot.cxx_
+  double mach = V_true / c;
   
   // publish Mach and TAS
   _mach_node->setDoubleValue(mach);
-  _tas_node->setDoubleValue(CS * mach);
+  _tas_node->setDoubleValue(V_true * SG_MPS_TO_KT );
 }
 
 // end of airspeed_indicator.cxx
diff --git a/src/Instrumentation/airspeed_indicator.hxx b/src/Instrumentation/airspeed_indicator.hxx
index 0a35406ae..1df6fe09b 100644
--- a/src/Instrumentation/airspeed_indicator.hxx
+++ b/src/Instrumentation/airspeed_indicator.hxx
@@ -1,6 +1,7 @@
 // airspeed_indicator.hxx - a regular VSI tied to the static port.
 // Written by David Megginson, started 2002.
 //
+// Last modified by Eric van den Berg, 24 Nov 2012
 // This file is in the Public Domain and comes with no warranty.
 
 
@@ -14,6 +15,9 @@
 #include <simgear/props/props.hxx>
 #include <simgear/structure/subsystem_mgr.hxx>
 
+// forward decls
+class FGEnvironmentMgr;
+
 /**
  * Model an airspeed indicator tied to the pitot and static ports.
  *
@@ -27,6 +31,8 @@
  * Output properties:
  *
  * /instrumentation/"name"/indicated-speed-kt
+ * /instrumentation/"name"/true-speed-kt
+ * /instrumentation/"name"/indicated-mach
  */
 class AirspeedIndicator : public SGSubsystem
 {
@@ -41,7 +47,7 @@ public:
     virtual void update (double dt);
 
 private:
-    void computeMach(double ias);
+    void computeMach();
 
     std::string _name;
     unsigned int _num;
@@ -65,8 +71,8 @@ private:
     SGPropertyNode_ptr _pressure_alt;
     SGPropertyNode_ptr _mach_node;
     SGPropertyNode_ptr _tas_node;
-    SGPropertyNode_ptr _sea_level_pressure_node;
-    SGPropertyNode_ptr _oat_celsius_node;
+    
+    FGEnvironmentMgr* _environmentManager;
 };
 
 #endif // __INSTRUMENTS_AIRSPEED_INDICATOR_HXX
diff --git a/src/Instrumentation/altimeter.cxx b/src/Instrumentation/altimeter.cxx
index 3ef204687..d9e858515 100644
--- a/src/Instrumentation/altimeter.cxx
+++ b/src/Instrumentation/altimeter.cxx
@@ -2,6 +2,7 @@
 // Written by David Megginson, started 2002.
 // Modified by John Denker in 2007 to use a two layer atmosphere
 // model in src/Environment/atmosphere.?xx
+// Last modified by Eric van den Berg, 25 Nov 2012
 //
 // This file is in the Public Domain and comes with no warranty.
 
@@ -19,6 +20,7 @@
 #  include <config.h>
 #endif
 
+#include <simgear/constants.h>
 #include <simgear/math/interpolater.hxx>
 #include <simgear/math/SGMath.hxx>
 
@@ -28,8 +30,6 @@
 
 #include "altimeter.hxx"
 
-const double hPa2inHg = 29.92 / 1013.25;
-
 Altimeter::Altimeter ( SGPropertyNode *node, double quantum )
     : _rootNode( 
        fgGetNode("/instrumentation",true)->
@@ -62,13 +62,13 @@ Altimeter::setSettingInHg( double value )
 double
 Altimeter::getSettingHPa() const
 {
-    return _settingInHg / hPa2inHg;
+    return _settingInHg * SG_INHG_TO_PA / 100;
 }
 
 void
 Altimeter::setSettingHPa( double value )
 {
-    _settingInHg = value * hPa2inHg;
+    _settingInHg = value * SG_PA_TO_INHG * 100;
 }
 
 
diff --git a/src/Systems/pitot.cxx b/src/Systems/pitot.cxx
index b367283b2..4f18431cd 100644
--- a/src/Systems/pitot.cxx
+++ b/src/Systems/pitot.cxx
@@ -1,6 +1,7 @@
 // pitot.cxx - the pitot air system.
 // Written by David Megginson, started 2002.
 //
+// Last modified by Eric van den Berg, 24 Nov 2012
 // This file is in the Public Domain and comes with no warranty.
 
 #ifdef HAVE_CONFIG_H
@@ -35,10 +36,9 @@ PitotSystem::init ()
     SGPropertyNode *node = fgGetNode(branch.c_str(), _num, true );
     _serviceable_node = node->getChild("serviceable", 0, true);
     _pressure_node = fgGetNode("/environment/pressure-inhg", true);
-    _density_node = fgGetNode("/environment/density-slugft3", true);
-    _velocity_node = fgGetNode("/velocities/airspeed-kt", true);
-    _slip_angle = fgGetNode("/orientation/side-slip-rad", true);
+    _mach_node = fgGetNode("/velocities/mach", true);
     _total_pressure_node = node->getChild("total-pressure-inhg", 0, true);
+    _measured_total_pressure_node = node->getChild("measured-total-pressure-inhg", 0, true);
 }
 
 void
@@ -51,29 +51,20 @@ PitotSystem::unbind ()
 {
 }
 
-#ifndef INHGTOPSF
-# define INHGTOPSF (2116.217/29.9212)
-#endif
-
-#ifndef PSFTOINHG
-# define PSFTOINHG (1/INHGTOPSF)
-#endif
-
-
 void
 PitotSystem::update (double dt)
 {
     if (_serviceable_node->getBoolValue()) {
-                                // The pitot tube sees the forward
-                                // velocity in the body axis.
-        double p = _pressure_node->getDoubleValue() * INHGTOPSF;
-        double r = _density_node->getDoubleValue();
-        double v = _velocity_node->getDoubleValue() * SG_KT_TO_FPS;
-        v *= cos(_slip_angle->getDoubleValue());
-        if (v < 0.0)
-            v = 0.0;
-        double q = 0.5 * r * v * v; // dynamic
-        _total_pressure_node->setDoubleValue((p + q) * PSFTOINHG);
+        double p = _pressure_node->getDoubleValue();
+        double mach = _mach_node->getDoubleValue();
+        mach = std::max( mach , 0.0 );
+        double p_t = p * pow(1 + 0.2 * mach*mach, 3.5 );    // true total pressure around aircraft
+        _total_pressure_node->setDoubleValue(p_t);
+        double p_t_meas = p_t;
+        if (mach > 1) {    
+          p_t_meas = p * pow( 1.2 * mach*mach, 3.5 ) * pow( 2.8/2.4*mach*mach - 0.4 / 2.4 , -2.5 );    // measured total pressure by pitot tube (Rayleigh formula, at Mach>1, normal shockwave in front of pitot tube)     
+        }
+        _measured_total_pressure_node->setDoubleValue(p_t_meas);
     }
 }
 
diff --git a/src/Systems/pitot.hxx b/src/Systems/pitot.hxx
index 1549c9eda..7f1e48787 100644
--- a/src/Systems/pitot.hxx
+++ b/src/Systems/pitot.hxx
@@ -1,6 +1,7 @@
 // pitot.hxx - the pitot air system.
 // Written by David Megginson, started 2002.
 //
+// Last modified by Eric van den Berg, 24 Nov 2012
 // This file is in the Public Domain and comes with no warranty.
 
 
@@ -29,13 +30,13 @@ using std::string;
  * Input properties:
  *
  * /systems/"name"/serviceable
- * /environment/pressure-slugft3
- * /environment/density-slugft3
- * /velocities/airspeed-kt
+ * /environment/pressure-inhg
+ * /velocities/mach
  *
  * Output properties:
  *
  * /systems/"name"/total-pressure-inhg
+ * /systems/"name"/measured-total-pressure-inhg
  */
 class PitotSystem : public SGSubsystem
 {
@@ -56,10 +57,9 @@ private:
     int _num;
     SGPropertyNode_ptr _serviceable_node;
     SGPropertyNode_ptr _pressure_node;
-    SGPropertyNode_ptr _density_node;
-    SGPropertyNode_ptr _velocity_node;
-    SGPropertyNode_ptr _slip_angle;
+    SGPropertyNode_ptr _mach_node;
     SGPropertyNode_ptr _total_pressure_node;
+    SGPropertyNode_ptr _measured_total_pressure_node;
 };
 
 #endif // __SYSTEMS_PITOT_HXX