1
0
Fork 0

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:
John Denker 2009-01-21 18:11:52 -07:00 committed by Tim Moore
parent 9a19c238a9
commit abc7ae96a5
5 changed files with 293 additions and 80 deletions

View file

@ -1,16 +1,49 @@
#include <simgear/math/SGMath.hxx>
#include <simgear/debug/logstream.hxx>
#include "atmosphere.hxx"
using namespace std;
#include <iostream>
FGAtmoCache::FGAtmoCache() :
a_tvs_p(0)
{}
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),
};
FGAtmoCache::~FGAtmoCache() {
delete a_tvs_p;
const int ISA_def_size(sizeof(ISA_def) / sizeof(ISA_layer));
#if 0
// Pressure within a layer, as a function of height
// and of layer parameters.
// Physics model: standard or nonstandard atmosphere,
// depending on what parameters you pass in.
// Heights in meters, pressures in Pa.
// As always, $lambda is positive in the troposphere,
// and zero in the first part of the stratosphere.
//
double Ph_layer(
const double hh,
const double hb,
const double Pb,
const double Tb,
const double lambda) {
using namespace atmodel;
if (lambda) {
// (1/N) &:=& {g\,m_M \over \lambda\,R}
double N = lambda * Rgas / mm / g;
return Pb * pow((Tb - lambda*(hh - hb)) / Tb, 1/N);
} else {
return Pb * exp(-g * mm / Rgas / Tb * (hh - hb));
}
}
// Pressure as a function of height.
@ -21,7 +54,7 @@ FGAtmoCache::~FGAtmoCache() {
//
// Height in meters, pressure in pascals.
double FGAtmo::p_vs_a(const double height) {
double old_p_vs_a(const double height) {
using namespace atmodel;
if (height <= 11000.) {
return P_layer(height, 0.0, ISA::P0, ISA::T0, ISA::lam0);
@ -33,10 +66,120 @@ double FGAtmo::p_vs_a(const double height) {
return 0;
}
// degrees C, height in feet
double FGAtmo::fake_t_vs_a_us(const double h_ft) {
#endif
// 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;
return ISA::T0 - ISA::lam0 * h_ft * foot - freezing;
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).
tuple<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_tuple(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_tuple(d0,d0);
}
FGAtmoCache::FGAtmoCache() :
a_tvs_p(0)
{}
FGAtmoCache::~FGAtmoCache() {
delete a_tvs_p;
}
/////////////
// 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,
const double Tsl) const {
using namespace atmodel;
return Tsl - ISA::lam0 * h_ft * foot;
}
// Dewpoint. degrees C or K, height in feet
@ -62,7 +205,8 @@ void FGAtmoCache::tabulate() {
a_tvs_p = new SGInterpTable;
for (double hgt = -1000; hgt <= 32000;) {
double press = p_vs_a(hgt);
double press,temp;
make_tuple(ref(press), ref(temp)) = PT_vs_hpt(hgt);
a_tvs_p->addEntry(press / inHg, hgt / foot);
#ifdef DEBUG_EXPORT_P_H
@ -86,25 +230,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 +257,8 @@ void FGAtmoCache::check_model() {
break;
using namespace atmodel;
cache();
double press = p_vs_a(height);
double press,temp;
make_tuple(ref(press), ref(temp)) = PT_vs_hpt(height);
cout << "Height: " << height
<< " \tpressure: " << press << endl;
cout << "Check: "
@ -154,21 +280,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 +323,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

@ -27,6 +27,8 @@
#include <simgear/compiler.h>
#include <simgear/math/interpolater.hxx>
#include "boost/tuple/tuple.hpp"
using namespace boost;
#include <cmath>
@ -50,8 +52,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 +69,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[];
tuple<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

@ -18,7 +18,7 @@
// along with this program; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
//
// $Id$
// $Id: environment.cxx,v 1.1 2009/01/30 15:07:04 jsd Exp $
#ifdef HAVE_CONFIG_H
@ -36,15 +36,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
@ -111,7 +114,7 @@ _setup_tables ()
atmosphere_data[i][2]);
}
}
#endif
////////////////////////////////////////////////////////////////////////
@ -137,7 +140,9 @@ void FGEnvironment::_init()
wind_from_down_fps = 0;
altitude_half_to_sun_m = 1000;
altitude_tropo_top_m = 10000;
#ifdef USING_TABLES
_setup_tables();
#endif
_recalc_density();
_recalc_relative_humidity();
@ -200,12 +205,12 @@ FGEnvironment::read (const SGPropertyNode * node)
&FGEnvironment::set_visibility_m);
if (!maybe_copy_value(this, node, "temperature-sea-level-degc",
&FGEnvironment::set_temperature_sea_level_degc))
&FGEnvironment::set_temperature_sea_level_degc))
maybe_copy_value(this, node, "temperature-degc",
&FGEnvironment::set_temperature_degc);
if (!maybe_copy_value(this, node, "dewpoint-sea-level-degc",
&FGEnvironment::set_dewpoint_sea_level_degc))
&FGEnvironment::set_dewpoint_sea_level_degc))
maybe_copy_value(this, node, "dewpoint-degc",
&FGEnvironment::set_dewpoint_degc);
@ -372,7 +377,7 @@ FGEnvironment::set_temperature_sea_level_degc (double t)
temperature_sea_level_degc = t;
if (dewpoint_sea_level_degc > t)
dewpoint_sea_level_degc = t;
_recalc_alt_temperature();
_recalc_alt_pt();
_recalc_density();
}
@ -381,6 +386,8 @@ FGEnvironment::set_temperature_degc (double t)
{
temperature_degc = t;
_recalc_sl_temperature();
_recalc_sl_pressure();
_recalc_alt_pt();
_recalc_density();
_recalc_relative_humidity();
}
@ -408,7 +415,7 @@ void
FGEnvironment::set_pressure_sea_level_inhg (double p)
{
pressure_sea_level_inhg = p;
_recalc_alt_pressure();
_recalc_alt_pt();
_recalc_density();
}
@ -471,9 +478,8 @@ void
FGEnvironment::set_elevation_ft (double e)
{
elevation_ft = e;
_recalc_alt_temperature();
_recalc_alt_dewpoint();
_recalc_alt_pressure();
_recalc_alt_pt();
_recalc_density();
_recalc_relative_humidity();
}
@ -538,26 +544,34 @@ FGEnvironment::_recalc_ne ()
sin(wind_from_heading_deg * SGD_DEGREES_TO_RADIANS);
}
// 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
@ -585,15 +599,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;
make_tuple(ref(press), ref(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

View file

@ -103,11 +103,11 @@ private:
void _recalc_ne ();
void _recalc_sl_temperature ();
void _recalc_alt_temperature ();
void _recalc_sl_dewpoint ();
void _recalc_alt_dewpoint ();
void _recalc_sl_pressure ();
void _recalc_alt_pressure ();
void _recalc_alt_pt ();
void _recalc_density ();
void _recalc_density_tropo_avg_kgm3 ();

View file

@ -718,6 +718,7 @@ struct GeneralInitOperation : public GraphicsContextOperation
static void fgIdleFunction ( void ) {
static osg::ref_ptr<GeneralInitOperation> genOp;
////////cerr << "Idle state: " << idle_state << endl;
if ( idle_state == 0 ) {
idle_state++;
// Pick some window on which to do queries.