1
0
Fork 0

Merge branch 'topic/atmos-merge' into next

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.
This commit is contained in:
Tim Moore 2009-09-12 15:50:06 +02:00
commit d11ad100cc
4 changed files with 259 additions and 93 deletions

View file

@ -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.
// 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);

View file

@ -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);
};

View file

@ -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
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();

View file

@ -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 ();