Merge branch 'jsd/atmos' into topic/atmos-merge

Conflicts:
	src/Environment/environment.cxx
	src/Environment/environment.hxx

John Denker's atmosphere changes. Original commit message:
Two-parameter physics-based model of atmosphere up to 262,467 ft i.e.
the top of the mesosphere. Correctly exhibits the HALT phenomenon.

src/Environment/atmosphere.cxx

@@ -1,10 +1,118 @@
+#include <boost/tuple/tuple.hpp>
+
 #include <simgear/math/SGMath.hxx>
+#include <simgear/debug/logstream.hxx>
 
 #include "atmosphere.hxx"
 
 using namespace std;
 #include <iostream>
 
+const ISA_layer ISA_def[] = {
+//        0    1        2        3           4         5      6         7         8
+//       id   (m)      (ft)     (Pa)        (inHg)    (K)     (C)      (K/m)     (K/ft)
+ISA_layer(0,      0,       0,   101325,    29.92126, 288.15,  15.00,  0.0065,   0.0019812),
+ISA_layer(1,  11000,   36089,  22632.1,    6.683246, 216.65, -56.50,       0,           0),
+ISA_layer(2,  20000,   65616,  5474.89,    1.616734, 216.65, -56.50, -0.0010,  -0.0003048),
+ISA_layer(3,  32000,  104986,  868.019,    0.256326, 228.65, -44.50, -0.0028,  -0.0008534),
+ISA_layer(4,  47000,  154199,  110.906,   0.0327506, 270.65,  -2.50,       0,           0),
+ISA_layer(5,  51000,  167322,  66.9389,   0.0197670, 270.65,  -2.50,  0.0028,   0.0008534),
+ISA_layer(6,  71000,  232939,  3.95642,  0.00116833, 214.65, -58.50,  0.0020,   0.0006096),
+ISA_layer(7,  80000,  262467,  0.88628, 0.000261718, 196.65, -76.50),
+};
+
+const int ISA_def_size(sizeof(ISA_def) / sizeof(ISA_layer));
+
+// Pressure within a layer, as a function of height.
+// Physics model:  standard or nonstandard atmosphere,
+//    depending on what parameters you pass in.
+// Height in meters, pressures in pascals.
+// As always, lapse is positive in the troposphere,
+// and zero in the first part of the stratosphere.
+
+double P_layer(const double height, const double href,
+        const double Pref, const double Tref,
+        const double lapse) {
+    using namespace atmodel;
+    if (lapse) {
+        double N = lapse * Rgas / mm / g;
+        return Pref * pow( (Tref - lapse*(height - href)) / Tref , (1/N));
+    } else {
+        return Pref * exp(-g * mm / Rgas / Tref * (height - href));
+    }
+}
+
+
+// Temperature within a layer, as a function of height.
+// Physics model:  standard or nonstandard atmosphere
+//  depending on what parameters you pass in.
+// $hh in meters, pressures in Pa.
+// As always, $lambda is positive in the troposphere,
+// and zero in the first part of the stratosphere.
+double T_layer (
+          const double hh, 
+          const double hb, 
+          const double Pb, 
+          const double Tb, 
+          const double lambda) {
+  return Tb - lambda*(hh - hb);
+}
+
+// Pressure and temperature as a function of height, Psl, and Tsl.
+// heights in meters, pressures in Pa.
+// Daisy chain version.
+// We need "seed" values for sea-level pressure and temperature.
+// In addition, for every layer, we need three things
+//  from the table: the reference height in that layer, 
+//  the lapse in that layer, and the cap (if any) for that layer 
+// (which we take from the /next/ row of the table, if any).
+pair<double,double> PT_vs_hpt(
+      const double hh, 
+      const double _p0,
+      const double _t0
+) {
+  
+  const double d0(0);
+  double hgt = ISA_def[0].height;
+  double p0 =  _p0;
+  double t0 =  _t0;
+#if 0
+    cout << "PT_vs_hpt: " << hh << "   " << p0 << "   " << t0 << endl;
+#endif 
+
+  int ii = 0;
+  for (const ISA_layer* pp = ISA_def; pp->lapse != -1; pp++, ii++) {
+#if 0
+    cout << "PT_vs_hpt: " << ii
+        << "  height: " << pp->height
+        << "  temp: "   << pp->temp
+        << "  lapse: "  << pp->lapse 
+        << endl;
+#endif
+    double xhgt(9e99);
+    double lapse = pp->lapse;
+// Stratosphere starts at a definite temperature,
+// not a definite height:
+    if (ii == 0) {
+      xhgt = hgt + (t0 - (pp+1)->temp) / lapse;
+    } else if ((pp+1)->lapse != -1) {
+      xhgt = (pp+1)->height;      
+    }
+    if (hh <= xhgt) {
+      return make_pair(P_layer(hh, hgt, p0, t0, lapse),
+                 T_layer(hh, hgt, p0, t0, lapse));
+    }
+    p0 = P_layer(xhgt, hgt, p0, t0, lapse);
+    t0 = t0 - lapse * (xhgt - hgt);
+    hgt = xhgt;
+  }
+  
+// Should never get here.
+  SG_LOG(SG_GENERAL, SG_ALERT, "PT_vs_hpt: ran out of layers");
+  return make_pair(d0, d0);
+}
+
+
 FGAtmoCache::FGAtmoCache() :
     a_tvs_p(0)
 {}
@@ -13,30 +121,19 @@ FGAtmoCache::~FGAtmoCache() {
     delete a_tvs_p;
 }
 
-// Pressure as a function of height.
-// Valid below 32000 meters,
-//   i.e. troposphere and first two layers of stratosphere.
-// Does not depend on any caching; can be used to
-// *construct* caches and interpolation tables.
-//
-// Height in meters, pressure in pascals.
-
-double FGAtmo::p_vs_a(const double height) {
-    using namespace atmodel;
-    if (height <= 11000.) {
-        return P_layer(height, 0.0, ISA::P0, ISA::T0, ISA::lam0);
-    } else if (height <= 20000.) {
-        return P_layer(height, 11000., 22632.06, 216.65, 0.0);
-    } else if (height <= 32000.) {
-        return P_layer(height, 20000.,  5474.89, 216.65, -0.001);
-    }
-    return 0;
-}
 
+/////////////
+// The following two routines are called "fake" because they
+// bypass the exceedingly complicated layer model implied by
+// the "weather conditioins" popup menu.
+// For now we must bypass it for several reasons, including
+// the fact that we don't have an "environment" object for
+// the airport (only for the airplane).
 // degrees C, height in feet
-double FGAtmo::fake_t_vs_a_us(const double h_ft) {
+double FGAtmo::fake_T_vs_a_us(const double h_ft, 
+                const double Tsl) const {
     using namespace atmodel;
-    return ISA::T0 - ISA::lam0 * h_ft * foot - freezing;
+    return Tsl - ISA::lam0 * h_ft * foot;
 }
 
 // Dewpoint.  degrees C or K, height in feet
@@ -62,7 +159,8 @@ void FGAtmoCache::tabulate() {
     a_tvs_p = new SGInterpTable;
 
     for (double hgt = -1000; hgt <= 32000;) {
-        double press = p_vs_a(hgt);
+        double press,temp;
+        boost::tie(press, temp) = PT_vs_hpt(hgt);
         a_tvs_p->addEntry(press / inHg, hgt / foot);
 
 #ifdef DEBUG_EXPORT_P_H
@@ -86,25 +184,6 @@ void FGAtmoCache::cache() {
         tabulate();
 }
 
-// Pressure within a layer, as a function of height.
-// Physics model:  standard or nonstandard atmosphere,
-//    depending on what parameters you pass in.
-// Height in meters, pressures in pascals.
-// As always, lapse is positive in the troposphere,
-// and zero in the first part of the stratosphere.
-
-double FGAtmo::P_layer(const double height, const double href,
-        const double Pref, const double Tref,
-        const double lapse) {
-    using namespace atmodel;
-    if (lapse) {
-        double N = lapse * Rgas / mm / g;
-        return Pref * pow( (Tref - lapse*(height - href)) / Tref , (1/N));
-    } else {
-        return Pref * exp(-g * mm / Rgas / Tref * (height - href));
-    }
-}
-
 // Check the basic function,
 // then compare against the interpolator.
 void FGAtmoCache::check_model() {
@@ -132,7 +211,8 @@ void FGAtmoCache::check_model() {
             break;
         using namespace atmodel;
         cache();
-        double press = p_vs_a(height);
+        double press,temp;
+        boost::tie(press, temp) = PT_vs_hpt(height);
         cout << "Height: " << height
              << " \tpressure: " << press << endl;
         cout << "Check:  "
@@ -154,21 +234,30 @@ double FGAltimeter::reading_ft(const double p_inHg, const double set_inHg) {
     return (press_alt - kollsman_shift);
 }
 
-// Altimeter setting.
+// Altimeter setting _in pascals_ 
-// Field elevation in feet
+//  ... caller gets to convert to inHg or millibars
-// Field pressure in inHg
+// Field elevation in m
-// field elevation in troposphere only
+// Field pressure in pascals
-double FGAtmo::qnh(const double field_ft, const double press_inHg) {
+// Valid for fields within the troposphere only.
+double FGAtmo::QNH(const double field_elev, const double field_press) {
     using namespace atmodel;
 
     //  Equation derived in altimetry.htm
     // exponent in QNH equation:
     double nn = ISA::lam0 * Rgas / g / mm;
     // pressure ratio factor:
-    double prat = pow(ISA::P0/inHg / press_inHg, nn);
+    double prat = pow(ISA::P0 / field_press, nn);
-
+    double rslt =  field_press
-    return press_inHg
+            * pow(1. + ISA::lam0 * field_elev / ISA::T0 * prat, 1./nn);
-            * pow(1 + ISA::lam0 * field_ft * foot / ISA::T0 * prat, 1/nn);
+#if 0
+    SG_LOG(SG_GENERAL, SG_ALERT, "QNH: elev: " << field_elev
+                << "  press: " << field_press
+                << "  prat: "  << prat
+                << "  rslt: " << rslt
+                << "  inHg: " << inHg
+                << "  rslt/inHG: " << rslt/inHg);
+#endif
+    return rslt;
 }
 
 void FGAltimeter::dump_stack1(const double Tref) {
@@ -188,13 +277,14 @@ void FGAltimeter::dump_stack1(const double Tref) {
     double hgts[] = {0, 2500, 5000, 7500, 10000, -9e99};
     for (int ii = 0; ; ii++) {
         double hgt_ft = hgts[ii];
+        double hgt = hgt_ft * foot;
         if (hgt_ft < -1e6)
             break;
-        double press = P_layer(hgt_ft*foot, 0, ISA::P0, Tref, ISA::lam0);
+        double press = P_layer(hgt, 0, ISA::P0, Tref, ISA::lam0);
-        double p_inHg = press / inHg;
+        double qnhx = QNH(hgt, press) / inHg;
-        double qnhx = qnh(hgt_ft, p_inHg);
         double qnh2 = SGMiscd::round(qnhx*100)/100;
 
+        double p_inHg = press / inHg;
         double Aprind = reading_ft(p_inHg);
         double Apr    = a_vs_p(p_inHg*inHg) / foot;
         double hind = reading_ft(p_inHg, qnh2);

src/Environment/atmosphere.hxx

@@ -29,6 +29,7 @@
 #include <simgear/math/interpolater.hxx>
 
 #include <cmath>
+#include <utility>
 
 using namespace std;
 
@@ -50,8 +51,9 @@ namespace atmodel {
     SCD(mm, .0289644);		// [kg/mole] molar mass of air (dry?)
     SCD(Rgas, 8.31432);		// [J/K/mole] gas constant
     SCD(inch, 0.0254);		// [m] definition of inch
-    SCD(foot, 12 * inch);		// [m]
+    SCD(foot, 12 * inch);	// [m]
     SCD(inHg, 101325.0 / 760 * 1000 * inch);  // [Pa] definition of inHg
+    SCD(mbar, 100.);            // [Pa] definition of millibar
     SCD(freezing, 273.15);	// [K] centigrade - kelvin offset
     SCD(nm, 1852);		// [m] nautical mile (NIST)
     SCD(sm, 5280*foot);		// [m] nautical mile (NIST)
@@ -66,19 +68,48 @@ namespace atmodel {
 
 
 
+class ISA_layer {
+public:
+  double height;
+  double temp;
+  double lapse;
+  ISA_layer(int, double h, double, double, double, double t, double,
+                    double l=-1, double=0)
+   :  height(h),        // [meters]
+      temp(t),          // [kelvin]
+      lapse(l)          // [K/m]
+  {}
+};
+
+extern const ISA_layer ISA_def[];
+
+std::pair<double,double> PT_vs_hpt(
+          const double hh, 
+          const double _p0 = atmodel::ISA::P0,
+          const double _t0 = atmodel::ISA::T0);
+
+double P_layer(const double height, const double href,
+  const double Pref, const double Tref, const double lapse );
+
+double T_layer(const double height, const double href,
+  const double Pref, const double Tref,  const double lapse );
+
 // The base class is little more than a namespace.
 // It has no constructor, no destructor, and no variables.
 class FGAtmo {
 public:
-    double p_vs_a(const double height);
     double a_vs_p(const double press, const double qnh = atmodel::ISA::P0);
-    double fake_t_vs_a_us(const double h_ft);
+    double fake_T_vs_a_us(const double h_ft, 
+                const double Tsl = atmodel::ISA::T0) const;
     double fake_dp_vs_a_us(const double dpsl, const double h_ft);
-    double P_layer(const double height, const double href,
-    const double Pref, const double Tref,
-    const double lapse );
     void check_one(const double height);
-    double qnh(const double field_ft, const double press_in);
+
+// Altimeter setting _in pascals_ 
+//  ... caller gets to convert to inHg or millibars
+// Field elevation in m
+// Field pressure in pascals
+// Valid for fields within the troposphere only.
+    double QNH(const double field_elev, const double field_press);
 };
 
 

src/Environment/environment.cxx

@@ -27,6 +27,10 @@
 
 #include <math.h>
 
+#include <boost/tuple/tuple.hpp>
+
+#include <plib/sg.h>
+
 #include <simgear/constants.h>
 #include <simgear/debug/logstream.hxx>
 #include <simgear/math/interpolater.hxx>
@@ -34,15 +38,18 @@
 #include <simgear/environment/visual_enviro.hxx>
 
 #include <Main/fg_props.hxx>
+#include <signal.h>
 
 #include "environment.hxx"
-
+#include "atmosphere.hxx"
 
 
 ////////////////////////////////////////////////////////////////////////
 // Atmosphere model.
 ////////////////////////////////////////////////////////////////////////
 
+#ifdef USING_TABLES
+
 // Calculated based on the ISA standard day, as found at e.g.
 // http://www.av8n.com/physics/altimetry.htm
 
@@ -109,7 +116,7 @@ _setup_tables ()
 				   atmosphere_data[i][2]);
   }
 }
-
+#endif
 
 
 ////////////////////////////////////////////////////////////////////////
@@ -137,7 +144,9 @@ void FGEnvironment::_init()
     ridge_lift_fps= 0;
     altitude_half_to_sun_m = 1000;
     altitude_tropo_top_m = 10000;
+#ifdef USING_TABLES
     _setup_tables();
+#endif
     _recalc_density();
     _recalc_relative_humidity();
     live_update = true;
@@ -233,9 +242,8 @@ FGEnvironment::read (const SGPropertyNode * node)
                      &FGEnvironment::set_turbulence_rate_hz);
     // calculate derived properties here to avoid duplicate expensive computations
     _recalc_ne();
-    _recalc_alt_temperature();
+    _recalc_alt_pt();
     _recalc_alt_dewpoint();
-    _recalc_alt_pressure();
     _recalc_density();
     _recalc_relative_humidity();
 
@@ -397,7 +405,7 @@ FGEnvironment::set_temperature_sea_level_degc (double t)
   if (dewpoint_sea_level_degc > t)
       dewpoint_sea_level_degc = t;
   if( live_update ) {
-    _recalc_alt_temperature();
+    _recalc_alt_pt();
     _recalc_density();
   }
 }
@@ -408,6 +416,8 @@ FGEnvironment::set_temperature_degc (double t)
   temperature_degc = t;
   if( live_update ) {
     _recalc_sl_temperature();
+    _recalc_sl_pressure();
+    _recalc_alt_pt();
     _recalc_density();
     _recalc_relative_humidity();
   }
@@ -441,7 +451,7 @@ FGEnvironment::set_pressure_sea_level_inhg (double p)
 {
   pressure_sea_level_inhg = p;
   if( live_update ) {
-    _recalc_alt_pressure();
+    _recalc_alt_pt();
     _recalc_density();
   }
 }
@@ -536,9 +546,8 @@ FGEnvironment::set_elevation_ft (double e)
 {
   elevation_ft = e;
   if( live_update ) {
-    _recalc_alt_temperature();
     _recalc_alt_dewpoint();
-    _recalc_alt_pressure();
+    _recalc_alt_pt();
     _recalc_density();
     _recalc_relative_humidity();
   }
@@ -614,26 +623,34 @@ FGEnvironment::_recalc_updraft ()
   wind_from_down_fps = thermal_lift_fps + ridge_lift_fps ;
 }
 
+// Intended to help with the interpretation of METAR data,
+// not for random in-flight outside-air temperatures.
 void
 FGEnvironment::_recalc_sl_temperature ()
 {
-  // If we're in the stratosphere, leave sea-level temp alone
-  if (elevation_ft < 38000) {
-    temperature_sea_level_degc = (temperature_degc + 273.15)
-        / _temperature_degc_table->interpolate(elevation_ft)
-      - 273.15;
-  }
-}
 
-void
+#if 0
-FGEnvironment::_recalc_alt_temperature ()
+  {
-{
+    SG_LOG(SG_GENERAL, SG_DEBUG, "recalc_sl_temperature: using "
-  if (elevation_ft < 38000) {
+      << temperature_degc << " @ " << elevation_ft << " :: " << this);
-    temperature_degc = (temperature_sea_level_degc + 273.15) *
-        _temperature_degc_table->interpolate(elevation_ft) - 273.15;
-  } else {
-    temperature_degc = -56.49;	// Stratosphere is constant
   }
+#endif
+
+  if (elevation_ft >= ISA_def[1].height) {
+    SG_LOG(SG_GENERAL, SG_ALERT, "recalc_sl_temperature: "
+        << "valid only in troposphere, not " << elevation_ft);
+    return;
+  }
+
+// Clamp: temperature of the stratosphere, in degrees C:
+  double t_strato = ISA_def[1].temp - atmodel::freezing;
+  if (temperature_degc < t_strato) temperature_sea_level_degc = t_strato;
+  else temperature_sea_level_degc = 
+      temperature_degc - elevation_ft * ISA_def[0].lapse;
+
+// Alternative implemenation:
+//  else temperature_sea_level_inhg = T_layer(0., elevation_ft * foot,
+//      pressure_inhg * inHg, temperature_degc + freezing, ISA_def[0].lapse) - freezing;
 }
 
 void
@@ -661,15 +678,45 @@ FGEnvironment::_recalc_alt_dewpoint ()
 void
 FGEnvironment::_recalc_sl_pressure ()
 {
-  pressure_sea_level_inhg =
+  using namespace atmodel;
-    pressure_inhg / _pressure_inhg_table->interpolate(elevation_ft);
+#if 0
+  {
+    SG_LOG(SG_GENERAL, SG_ALERT, "recalc_sl_pressure: using "
+      << pressure_inhg << " and "
+      << temperature_degc << " @ " << elevation_ft << " :: " << this);
+  }
+#endif
+  pressure_sea_level_inhg = P_layer(0., elevation_ft * foot,
+      pressure_inhg * inHg, temperature_degc + freezing, ISA_def[0].lapse) / inHg;
 }
 
+// This gets called at frame rate, to account for the aircraft's
+// changing altitude. 
+// Called by set_elevation_ft() which is called by FGEnvironmentMgr::update
+
 void
-FGEnvironment::_recalc_alt_pressure ()
+FGEnvironment::_recalc_alt_pt ()
 {
-  pressure_inhg =
+  using namespace atmodel;
-    pressure_sea_level_inhg * _pressure_inhg_table->interpolate(elevation_ft);
+#if 0
+  {
+    static int count(0);
+    if (++count % 1000 == 0) {
+      SG_LOG(SG_GENERAL, SG_ALERT, 
+           "recalc_alt_pt for: " << elevation_ft
+        << "  using "  << pressure_sea_level_inhg 
+        << "  and "  << temperature_sea_level_degc
+        << " :: " << this
+        << "  # " << count);
+        ///////////////////////////////////raise(SIGUSR1);
+    }
+  }
+#endif
+  double press, temp;
+  boost::tie(press, temp) = PT_vs_hpt(elevation_ft * foot, 
+        pressure_sea_level_inhg * inHg, temperature_sea_level_degc + freezing);
+  temperature_degc = temp - freezing;
+  pressure_inhg = press / inHg;
 }
 
 void
@@ -814,9 +861,8 @@ interpolate (const FGEnvironment * env1, const FGEnvironment * env2,
 
     // calculate derived properties here to avoid duplicate expensive computations
     result->_recalc_ne();
-    result->_recalc_alt_temperature();
+    result->_recalc_alt_pt();
     result->_recalc_alt_dewpoint();
-    result->_recalc_alt_pressure();
     result->_recalc_density();
     result->_recalc_relative_humidity();
 

src/Environment/environment.hxx

@@ -104,11 +104,10 @@ public:
   virtual bool set_live_update(bool live_update);
 
   void _recalc_ne ();
-  void _recalc_alt_temperature ();
   void _recalc_alt_dewpoint ();
-  void _recalc_alt_pressure ();
   void _recalc_density ();
   void _recalc_relative_humidity ();
+  void _recalc_alt_pt ();
 private:
   void _init();
   void _recalc_hdgspd ();