diff --git a/src/Environment/atmosphere.cxx b/src/Environment/atmosphere.cxx index a15c943ed..099669586 100644 --- a/src/Environment/atmosphere.cxx +++ b/src/Environment/atmosphere.cxx @@ -260,6 +260,21 @@ double FGAtmo::QNH(const double field_elev, const double field_press) { return rslt; } +// 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) { using namespace atmodel; const int bs(200); diff --git a/src/Environment/atmosphere.hxx b/src/Environment/atmosphere.hxx index 55ccfb3f9..1ee4e1b6b 100644 --- a/src/Environment/atmosphere.hxx +++ b/src/Environment/atmosphere.hxx @@ -110,6 +110,15 @@ public: // Field pressure in pascals // Valid for fields within the troposphere only. double QNH(const double field_elev, const double field_press); +/** + * Invert the QNH calculation to get the field pressure from a metar + * report. Valid for fields within the troposphere only. + * @param field_elev field elevation in m + * @param qnh altimeter setting in pascals + * @return field pressure _in pascals_. Caller gets to convert to inHg + * or millibars + */ + static double fieldPressure(const double field_elev, const double qnh); }; diff --git a/src/Environment/environment_ctrl.cxx b/src/Environment/environment_ctrl.cxx index 5b9cfcfb9..330f8e010 100644 --- a/src/Environment/environment_ctrl.cxx +++ b/src/Environment/environment_ctrl.cxx @@ -34,6 +34,7 @@ #include <Main/fg_props.hxx> #include <Main/util.hxx> +#include "atmosphere.hxx" #include "fgmetar.hxx" #include "environment_ctrl.hxx" @@ -433,6 +434,25 @@ static inline double convert_to_180( double d ) return d > 180.0 ? d - 360.0 : d; } +// Return the sea level pressure for a metar observation, in inHg. +// This is different from QNH because it accounts for the current +// temperature at the observation point. +// metarPressure in inHg +// fieldHt in ft +// fieldTemp in C + +static double reducePressureSl(double metarPressure, double fieldHt, + double fieldTemp) +{ + double elev = fieldHt * SG_FEET_TO_METER; + double fieldPressure + = FGAtmo::fieldPressure(elev, metarPressure * atmodel::inHg); + double slPressure = P_layer(0, elev, fieldPressure, + fieldTemp + atmodel::freezing, + atmodel::ISA::lam0); + return slPressure / atmodel::inHg; +} + void FGMetarCtrl::update(double dt) { @@ -444,6 +464,7 @@ FGMetarCtrl::update(double dt) bool reinit_required = false; bool layer_rebuild_required = false; + double station_elevation_ft = station_elevation_n->getDoubleValue(); if (first_update) { double dir = base_wind_dir_n->getDoubleValue()+magnetic_variation_n->getDoubleValue(); @@ -455,7 +476,10 @@ FGMetarCtrl::update(double dt) fgDefaultWeatherValue("visibility-m", metarvis); double metarpressure = pressure_n->getDoubleValue(); - fgDefaultWeatherValue("pressure-sea-level-inhg", metarpressure); + fgDefaultWeatherValue("pressure-sea-level-inhg", + reducePressureSl(metarpressure, + station_elevation_ft, + temperature_n->getDoubleValue())); // We haven't already loaded a METAR, so apply it immediately. vector<SGPropertyNode_ptr> layers = clouds_n->getChildren("layer"); @@ -597,8 +621,11 @@ FGMetarCtrl::update(double dt) double pressure = boundary_sea_level_pressure_n->getDoubleValue(); double metarpressure = pressure_n->getDoubleValue(); - if( pressure != metarpressure ) { - pressure = interpolate_val( pressure, metarpressure, MaxPressureChangeInHgSec ); + double newpressure = reducePressureSl(metarpressure, + station_elevation_ft, + temperature_n->getDoubleValue()); + if( pressure != newpressure ) { + pressure = interpolate_val( pressure, newpressure, MaxPressureChangeInHgSec ); fgDefaultWeatherValue("pressure-sea-level-inhg", pressure); reinit_required = true; } @@ -666,11 +693,8 @@ FGMetarCtrl::update(double dt) } } } - { - double station_elevation_ft = station_elevation_n->getDoubleValue(); - set_temp_at_altitude(temperature_n->getDoubleValue(), station_elevation_ft); - set_dewpoint_at_altitude(dewpoint_n->getDoubleValue(), station_elevation_ft); - } + set_temp_at_altitude(temperature_n->getDoubleValue(), station_elevation_ft); + set_dewpoint_at_altitude(dewpoint_n->getDoubleValue(), station_elevation_ft); //TODO: check if temperature/dewpoint have changed. This requires reinit. // Force an update of the 3D clouds