X-Git-Url: https://git.mxchange.org/?a=blobdiff_plain;f=src%2FEnvironment%2Fatmosphere.cxx;h=b7130aea3e9a4a65cce7c4831cfaa49df985972c;hb=1edffd1d10bbe30a418e88140df54791ced23649;hp=2e41c0570e8dc09c21761f9339ce88a7c085721d;hpb=724ba5b4ad5cde75e49195fa6ef6f6cb642b3b23;p=flightgear.git diff --git a/src/Environment/atmosphere.cxx b/src/Environment/atmosphere.cxx index 2e41c0570..b7130aea3 100644 --- a/src/Environment/atmosphere.cxx +++ b/src/Environment/atmosphere.cxx @@ -1,9 +1,118 @@ +#include + #include +#include #include "atmosphere.hxx" using namespace std; #include +#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_ENVIRONMENT, SG_ALERT, "PT_vs_hpt: ran out of layers for h=" << hh ); + return make_pair(d0, d0); +} + FGAtmoCache::FGAtmoCache() : a_tvs_p(0) @@ -13,30 +122,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 +160,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 +185,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 +212,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: " @@ -142,7 +223,7 @@ void FGAtmoCache::check_model() { ////////////////////////////////////////////////////////////////////// -FGAltimeter::FGAltimeter() : kset(atmodel::ISA::P0), kft(0) +FGAltimeter::FGAltimeter() { cache(); } @@ -154,21 +235,45 @@ 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); + 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_ENVIRONMENT, SG_ALERT, "QNH: elev: " << field_elev + << " press: " << field_press + << " prat: " << prat + << " rslt: " << rslt + << " inHg: " << inHg + << " rslt/inHG: " << rslt/inHg); +#endif + return rslt; +} - return press_inHg - * pow(1 + ISA::lam0 * field_ft * foot / ISA::T0 * prat, 1/nn); +// Invert the QNH calculation to get the field pressure from a metar +// report. +// field pressure _in pascals_ +// ... caller gets to convert to inHg or millibars +// Field elevation in m +// Altimeter setting (QNH) in pascals +// Valid for fields within the troposphere only. +double FGAtmo::fieldPressure(const double field_elev, const double qnh) +{ + using namespace atmodel; + static const double nn = ISA::lam0 * Rgas / g / mm; + const double pratio = pow(qnh / ISA::P0, nn); + return ISA::P0 * pow(pratio - field_elev * ISA::lam0 / ISA::T0, 1.0 / nn); } void FGAltimeter::dump_stack1(const double Tref) { @@ -188,13 +293,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);