#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.
//
// 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);
return 0;
}
+
+#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;
+ 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) {
+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
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
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() {
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: "
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) {
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);
#include <simgear/compiler.h>
#include <simgear/math/interpolater.hxx>
+#include "boost/tuple/tuple.hpp"
+using namespace boost;
#include <cmath>
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)
+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);
};
// 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
#include <simgear/environment/visual_enviro.hxx>
#include <Main/fg_props.hxx>
+#include <signal.h>
#include "environment.hxx"
-
+#include "atmosphere.hxx"
\f
////////////////////////////////////////////////////////////////////////
// Atmosphere model.
////////////////////////////////////////////////////////////////////////
+#ifdef USING_TABLES
+
// Calculated based on the ISA standard day, as found at e.g.
// http://www.av8n.com/physics/altimetry.htm
atmosphere_data[i][2]);
}
}
-
+#endif
\f
////////////////////////////////////////////////////////////////////////
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();
&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);
temperature_sea_level_degc = t;
if (dewpoint_sea_level_degc > t)
dewpoint_sea_level_degc = t;
- _recalc_alt_temperature();
+ _recalc_alt_pt();
_recalc_density();
}
{
temperature_degc = t;
_recalc_sl_temperature();
+ _recalc_sl_pressure();
+ _recalc_alt_pt();
_recalc_density();
_recalc_relative_humidity();
}
FGEnvironment::set_pressure_sea_level_inhg (double p)
{
pressure_sea_level_inhg = p;
- _recalc_alt_pressure();
+ _recalc_alt_pt();
_recalc_density();
}
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();
}
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;
+
+#if 0
+ {
+ SG_LOG(SG_GENERAL, SG_DEBUG, "recalc_sl_temperature: using "
+ << temperature_degc << " @ " << elevation_ft << " :: " << this);
}
-}
+#endif
-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 (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
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
projects / flightgear.git / commitdiff