diff --git a/src/Environment/atmosphere.cxx b/src/Environment/atmosphere.cxx index 2e41c0570..a15c943ed 100644 --- a/src/Environment/atmosphere.cxx +++ b/src/Environment/atmosphere.cxx @@ -1,10 +1,118 @@ +#include + #include +#include #include "atmosphere.hxx" using namespace std; #include +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 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. -// Field elevation in feet -// Field pressure in inHg -// field elevation in troposphere only -double FGAtmo::qnh(const double field_ft, const double press_inHg) { +// 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 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); - - return press_inHg - * pow(1 + ISA::lam0 * field_ft * foot / ISA::T0 * prat, 1/nn); + double prat = pow(ISA::P0 / field_press, nn); + double rslt = field_press + * pow(1. + ISA::lam0 * field_elev / 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 p_inHg = press / inHg; - double qnhx = qnh(hgt_ft, p_inHg); + double press = P_layer(hgt, 0, ISA::P0, Tref, ISA::lam0); + double qnhx = QNH(hgt, press) / 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); diff --git a/src/Environment/atmosphere.hxx b/src/Environment/atmosphere.hxx index 12df63327..55ccfb3f9 100644 --- a/src/Environment/atmosphere.hxx +++ b/src/Environment/atmosphere.hxx @@ -29,6 +29,7 @@ #include #include +#include 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 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); }; diff --git a/src/Environment/environment.cxx b/src/Environment/environment.cxx index 74e1509dc..f915e5261 100644 --- a/src/Environment/environment.cxx +++ b/src/Environment/environment.cxx @@ -27,6 +27,10 @@ #include +#include + +#include + #include #include #include @@ -34,15 +38,18 @@ #include #include
+#include #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 -FGEnvironment::_recalc_alt_temperature () -{ - if (elevation_ft < 38000) { - temperature_degc = (temperature_sea_level_degc + 273.15) * - _temperature_degc_table->interpolate(elevation_ft) - 273.15; - } else { - temperature_degc = -56.49; // Stratosphere is constant +#if 0 + { + SG_LOG(SG_GENERAL, SG_DEBUG, "recalc_sl_temperature: using " + << temperature_degc << " @ " << elevation_ft << " :: " << this); } +#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 = - pressure_inhg / _pressure_inhg_table->interpolate(elevation_ft); + using namespace atmodel; +#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 = - pressure_sea_level_inhg * _pressure_inhg_table->interpolate(elevation_ft); + using namespace atmodel; +#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(); diff --git a/src/Environment/environment.hxx b/src/Environment/environment.hxx index 7e3e5b4b2..d9543d5ab 100644 --- a/src/Environment/environment.hxx +++ b/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 ();