Generate sea level pressure from metar

This needs to account for the current temperature.

diff --git a/src/Environment/atmosphere.cxx b/src/Environment/atmosphere.cxx
index a15c943ed47aba514bddb6212459f38f9a85c3fc..099669586a82b41f046ef5a7ecedb91d58be26ca 100644 (file)
--- 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 55ccfb3f9e580a76d6d08c4ee1006fcb714bef02..1ee4e1b6beeaa58da610de8d854737f0b31f4bc5 100644 (file)
--- 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 5b9cfcfb99744b5e0a72843dbb252f57ab1518d8..330f8e010152a57cfd2b66e0014bf77fab5e1d6b 100644 (file)
--- 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