Move aircraft-performance code into a public API

Upcoming flight-planning changes want to use the perf computation code
in route-path, so move it to a public class, and implement some of the
missing functionality, especially correct GS computation for altitude.

src/Aircraft/AircraftPerformance.cxx

@@ -0,0 +1,390 @@
+// AircraftPerformance.cxx - compute data about planned acft performance
+//
+// Copyright (C) 2018  James Turner  <james@flightgear.org>
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of the
+// License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+
+#include "AircraftPerformance.hxx"
+
+#include <cassert>
+
+#include <simgear/constants.h>
+
+#include <Main/fg_props.hxx>
+
+using namespace flightgear;
+
+double distanceForTimeAndSpeeds(double tSec, double v1, double v2)
+{
+    return tSec * 0.5 * (v1 + v2);
+}
+
+AircraftPerformance::AircraftPerformance()
+{
+    // read aircraft supplied performance data
+    if (fgGetNode("/aircraft/performance/bracket")) {
+        readPerformanceData();
+    } else {
+        // falls back to heuristic determination of the category,
+        // and a plausible default
+        icaoCategoryData();
+    }
+}
+
+double AircraftPerformance::groundSpeedForAltitudeKnots(int altitudeFt) const
+{
+    auto bracket = bracketForAltitude(altitudeFt);
+    return bracket->gsForAltitude(altitudeFt);
+}
+
+int AircraftPerformance::computePreviousAltitude(double distanceM, int targetAltFt) const
+{
+    auto bracket = bracketForAltitude(targetAltFt);
+    auto d = bracket->descendDistanceM(bracket->atOrBelowAltitudeFt, targetAltFt);
+    if (d < distanceM) {
+        // recurse to previous bracket
+        return computePreviousAltitude(distanceM - d, bracket->atOrBelowAltitudeFt+1);
+    }
+
+    // work out how far we travel laterally per foot change in altitude
+    // this value is in metres, we have to map FPM and GS in Knots to make
+    // everything work out
+    const double gsMPS =  bracket->gsForAltitude(targetAltFt) * SG_KT_TO_MPS;
+    const double t = distanceM / gsMPS;
+    return targetAltFt + bracket->descentRateFPM * (t / 60.0);
+}
+
+int AircraftPerformance::computeNextAltitude(double distanceM, int initialAltFt) const
+{
+    auto bracket = bracketForAltitude(initialAltFt);
+    auto d = bracket->climbDistanceM(initialAltFt, bracket->atOrBelowAltitudeFt);
+    if (d < distanceM) {
+        // recurse to next bracket
+        return computeNextAltitude(distanceM - d, bracket->atOrBelowAltitudeFt+1);
+    }
+
+    // work out how far we travel laterally per foot change in altitude
+    // this value is in metres, we have to map FPM and GS in Knots to make
+    // everything work out
+    const double gsMPS =  bracket->gsForAltitude(initialAltFt) * SG_KT_TO_MPS;
+    const double t = distanceM / gsMPS;
+    return initialAltFt + bracket->climbRateFPM * (t / 60.0);
+}
+
+static string_list readTags()
+{
+    string_list r;
+    const auto tagsNode = fgGetNode("/sim/tags");
+    if (!tagsNode)
+        return r;
+    
+    for (auto t : tagsNode->getChildren("tag")) {
+        r.push_back(t->getStringValue());
+    }
+    return r;
+}
+
+static bool stringListContains(const string_list& t, const std::string& s)
+{
+    auto it = std::find(t.begin(), t.end(), s);
+    return it != t.end();
+}
+
+std::string AircraftPerformance::heuristicCatergoryFromTags() const
+{
+    const auto tags(readTags());
+    
+    if (stringListContains(tags, "turboprop"))
+        return {ICAO_AIRCRAFT_CATEGORY_C};
+    
+    // any way we could distuinguish fast and slow GA aircraft?
+    if (stringListContains(tags, "ga")) {
+        return {ICAO_AIRCRAFT_CATEGORY_A};
+    }
+    
+    if (stringListContains(tags, "jet")) {
+        return {ICAO_AIRCRAFT_CATEGORY_E};
+    }
+    
+    return {ICAO_AIRCRAFT_CATEGORY_C};
+}
+
+void AircraftPerformance::icaoCategoryData()
+{
+    std::string propCat = fgGetString("/aircraft/performance/icao-category");
+    if (propCat.empty()) {
+        propCat = heuristicCatergoryFromTags();
+    }
+
+    const char aircraftCategory = propCat.front();
+     //     pathTurnRate = 3.0; // 3 deg/sec = 180deg/min = standard rate turn
+      switch (aircraftCategory) {
+      case ICAO_AIRCRAFT_CATEGORY_A:
+          _perfData.push_back(Bracket(4000, 600, 1200, 75));
+          _perfData.push_back(Bracket(10000, 600, 1200, 140));
+          break;
+
+      case ICAO_AIRCRAFT_CATEGORY_B:
+          _perfData.push_back(Bracket(4000, 100, 1200, 100));
+          _perfData.push_back(Bracket(10000, 800, 1200, 160));
+          _perfData.push_back(Bracket(18000, 600, 1800, 200));
+          break;
+
+      case ICAO_AIRCRAFT_CATEGORY_C:
+          _perfData.push_back(Bracket(4000, 1800, 1800, 150));
+          _perfData.push_back(Bracket(10000, 1800, 1800, 200));
+          _perfData.push_back(Bracket(18000, 1200, 1800, 270));
+          _perfData.push_back(Bracket(60000, 800, 1200, 0.80, true /* is Mach */));
+          break;
+
+      case ICAO_AIRCRAFT_CATEGORY_D:
+      case ICAO_AIRCRAFT_CATEGORY_E:
+      default:
+          _perfData.push_back(Bracket(4000, 1800, 1800, 180));
+          _perfData.push_back(Bracket(10000, 1800, 1800, 230));
+          _perfData.push_back(Bracket(18000, 1200, 1800, 270));
+          _perfData.push_back(Bracket(60000, 800, 1200, 0.87, true /* is Mach */));
+          break;
+      }
+}
+
+void AircraftPerformance::readPerformanceData()
+{
+    for (auto nd : fgGetNode("/aircraft/performance/")->getChildren("bracket")) {
+        const int atOrBelowAlt = nd->getIntValue("at-or-below-ft");
+        const int climbFPM = nd->getIntValue("climb-rate-fpm");
+        const int descentFPM = nd->getIntValue("descent-rate-fpm");
+        bool isMach = nd->hasChild("speed-mach");
+        double speed;
+        if (isMach) {
+            speed = nd->getDoubleValue("speed-mach");
+        } else {
+            speed = nd->getIntValue("speed-ias-knots");
+        }
+
+        Bracket b(atOrBelowAlt, climbFPM, descentFPM, speed, isMach);
+        _perfData.push_back(b);
+    }
+}
+
+auto AircraftPerformance::bracketForAltitude(int altitude) const
+        -> PerformanceVec::const_iterator
+{
+    assert(!_perfData.empty());
+    if (_perfData.front().atOrBelowAltitudeFt >= altitude)
+        return _perfData.begin();
+
+    for (auto it = _perfData.begin(); it != _perfData.end(); ++it) {
+        if (it->atOrBelowAltitudeFt > altitude) {
+            return it;
+        }
+    }
+
+    return _perfData.end() - 1;
+}
+
+auto AircraftPerformance::rangeForAltitude(int lowAltitude, int highAltitude) const
+        -> BracketRange
+{
+    return {bracketForAltitude(lowAltitude), bracketForAltitude(highAltitude)};
+}
+
+void AircraftPerformance::traverseAltitudeRange(int initialElevationFt, int targetElevationFt,
+                                                TraversalFunc tf) const
+{
+    auto r = rangeForAltitude(initialElevationFt, targetElevationFt);
+    if (r.first == r.second) {
+       tf(*r.first, initialElevationFt, targetElevationFt);
+       return;
+    }
+
+    if (initialElevationFt < targetElevationFt) {
+        tf(*r.first, initialElevationFt, r.first->atOrBelowAltitudeFt);
+        int previousBracketCapAltitude = r.first->atOrBelowAltitudeFt;
+        for (auto bracket = r.first + 1; bracket != r.second; ++bracket) {
+            tf(*bracket, previousBracketCapAltitude, bracket->atOrBelowAltitudeFt);
+            previousBracketCapAltitude = bracket->atOrBelowAltitudeFt;
+        }
+        
+        tf(*r.second, previousBracketCapAltitude, targetElevationFt);
+    } else {
+        int nextBracketCapAlt = (r.first - 1)->atOrBelowAltitudeFt;
+        tf(*r.first, initialElevationFt, nextBracketCapAlt);
+        for (auto bracket = r.first - 1; bracket != r.second; --bracket) {
+            nextBracketCapAlt = (r.first - 1)->atOrBelowAltitudeFt;
+            tf(*bracket, bracket->atOrBelowAltitudeFt, nextBracketCapAlt);
+        }
+        
+        tf(*r.second, nextBracketCapAlt, targetElevationFt);
+    }
+}
+
+double AircraftPerformance::distanceNmBetween(int initialElevationFt, int targetElevationFt) const
+{
+    double result = 0.0;
+    TraversalFunc tf = [&result](const Bracket& bk, int alt1, int alt2) {
+        result += (alt1 > alt2) ? bk.descendDistanceM(alt1, alt2) : bk.climbDistanceM(alt1, alt2);
+    };
+    traverseAltitudeRange(initialElevationFt, targetElevationFt, tf);
+    return result * SG_METER_TO_NM;
+}
+
+double AircraftPerformance::timeBetween(int initialElevationFt, int targetElevationFt) const
+{
+    double result = 0.0;
+    TraversalFunc tf = [&result](const Bracket& bk, int alt1, int alt2) {
+        SG_LOG(SG_GENERAL, SG_INFO, "Range:" << alt1 << " " << alt2);
+        result += (alt1 > alt2) ? bk.descendTime(alt1, alt2) : bk.climbTime(alt1, alt2);
+    };
+    traverseAltitudeRange(initialElevationFt, targetElevationFt, tf);
+    return result;
+}
+
+double AircraftPerformance::timeToCruise(double cruiseDistanceNm, int cruiseAltitudeFt) const
+{
+    auto b = bracketForAltitude(cruiseAltitudeFt);
+    return (cruiseDistanceNm / b->gsForAltitude(cruiseAltitudeFt)) * 3600.0;
+}
+
+double oatCForAltitudeFt(int altitudeFt)
+{
+    if (altitudeFt > 36089)
+        return -56.5;
+
+    // lapse rate in C per ft
+    const double T_r = .0019812;
+    return 15.0 - (altitudeFt * T_r);
+}
+
+double oatKForAltitudeFt(int altitudeFt)
+{
+    return oatCForAltitudeFt(altitudeFt) + 273.15;
+}
+
+double pressureAtAltitude(int altitude)
+{
+#if 0
+  p= P_0*(1-6.8755856*10^-6 h)^5.2558797    h<36,089.24ft
+  p_Tr= 0.2233609*P_0                  
+  p=p_Tr*exp(-4.806346*10^-5(h-36089.24)) h>36,089.24ft 
+
+    magic numbers
+    6.8755856*10^-6 = T'/T_0, where T' is the standard temperature lapse rate and T_0 is the standard sea-level temperature.
+    5.2558797 = Mg/RT', where M is the (average) molecular weight of air, g is the acceleration of gravity and R is the gas constant.
+    4.806346*10^-5 = Mg/RT_tr, where T_tr is the temperature at the tropopause.
+#endif
+
+    const double k = 6.8755856e-6;
+    const double MgRT = 5.2558797;
+    const double MgRT_tr = 4.806346e-5;
+    const double P_0 = 29.92126; // (standard) sea-level pressure
+    if (altitude > 36089) {
+        const double P_Tr = 0.2233609 * P_0;
+        const double altAboveTr = altitude - 36089;
+        return P_Tr * exp(MgRT_tr * altAboveTr);
+    } else {
+        return P_0 * pow(1.0 - (k * altitude), MgRT);
+    }
+}
+
+double computeMachFromIAS(int iasKnots, int altitudeFt)
+{
+#if 0
+ // from the aviation formulary
+  DP=P_0*((1 + 0.2*(IAS/CS_0)^2)^3.5 -1)
+  M=(5*( (DP/P + 1)^(2/7) -1) )^0.5   (*)
+#endif
+    const double Cs_0 = 661.4786; // speed of sound at sea level, knots
+    const double P_0 = 29.92126; // (standard) sea-level pressure
+    const double iasCsRatio = iasKnots / Cs_0;
+    const double P = pressureAtAltitude(altitudeFt);
+    // differential pressure
+    const double DP = P_0 * (pow(1.0 + 0.2 * pow(iasCsRatio, 2.0), 3.5) - 1.0);
+
+    const double pressureRatio = DP / P + 1.0;
+    const double M = pow(5.0 * (pow(pressureRatio, 2.0 / 7.0) - 1.0), 0.5);
+    if (M > 1.0) {
+        SG_LOG(SG_GENERAL, SG_INFO, "computeMachFromIAS: computed Mach is supersonic, fix for shock wave");
+    }
+    return M;
+}
+
+int AircraftPerformance::Bracket::gsForAltitude(int altitude) const
+{
+    double M = 0.0;
+    if (speedIsMach) {
+        M = speedIASOrMach; // simple
+    } else {
+        M = computeMachFromIAS(speedIASOrMach, altitude);
+    }
+
+    // CS = sound speed= 38.967854*sqrt(T+273.15)  where T is the OAT in celsius.
+    const double CS = 38.967854 * sqrt(oatKForAltitudeFt(altitude));
+    const double TAS = M * CS;
+    return TAS;
+}
+
+double AircraftPerformance::Bracket::climbTime(int alt1, int alt2) const
+{
+    return (alt2 - alt1) / static_cast<double>(climbRateFPM) * 60.0;
+}
+
+double AircraftPerformance::Bracket::climbDistanceM(int alt1, int alt2) const
+{
+    const double t = climbTime(alt1, alt2);
+    return distanceForTimeAndSpeeds(t,
+                                    SG_KT_TO_MPS * gsForAltitude(alt1),
+                                    SG_KT_TO_MPS * gsForAltitude(alt2));
+}
+
+double AircraftPerformance::Bracket::descendTime(int alt1, int alt2) const
+{
+    return (alt1 - alt2) / static_cast<double>(descentRateFPM) * 60.0;
+}
+
+double AircraftPerformance::Bracket::descendDistanceM(int alt1, int alt2) const
+{
+    const double t = descendTime(alt1, alt2);
+    return distanceForTimeAndSpeeds(t,
+                                    SG_KT_TO_MPS * gsForAltitude(alt1),
+                                    SG_KT_TO_MPS * gsForAltitude(alt2));
+}
+
+double AircraftPerformance::turnRadiusMForAltitude(int altitudeFt) const
+{
+#if 0
+    From the aviation formulary again
+    In a steady turn, in no wind, with bank angle, b at an airspeed v
+    
+    tan(b)= v^2/(R g)
+    
+    With R in feet, v in knots, b in degrees and w in degrees/sec (inconsistent units!), numerical constants are introduced:
+    
+    R =v^2/(11.23*tan(0.01745*b))
+    (Example) At 100 knots, with a 45 degree bank, the radius of turn is 100^2/(11.23*tan(0.01745*45))= 891 feet.
+    
+    The bank angle b_s for a standard rate turn is given by:
+        
+        b_s = 57.3*atan(v/362.1)
+        (Example) for 100 knots, b_s = 57.3*atan(100/362.1) = 15.4 degrees
+            
+    Working in meter-per-second and radians removes a bunch of constants again.
+#endif
+    const double gsKts = groundSpeedForAltitudeKnots(altitudeFt);
+    const double gs = gsKts * SG_KT_TO_MPS;
+    const double bankAngleRad = atan(gsKts/362.1);
+    const double r = (gs * gs)/(SG_g0_m_p_s2 * tan(bankAngleRad));
+    return r;
+}
+

src/Aircraft/AircraftPerformance.hxx

@@ -0,0 +1,112 @@
+// AircraftPerformance.hxx - compute data about planned acft performance
+//
+// Copyright (C) 2018  James Turner  <james@flightgear.org>
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of the
+// License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+
+#ifndef AIRCRAFTPERFORMANCE_HXX
+#define AIRCRAFTPERFORMANCE_HXX
+
+#include <vector>
+#include <functional>
+
+namespace flightgear
+{
+
+const char ICAO_AIRCRAFT_CATEGORY_A = 'A';
+const char ICAO_AIRCRAFT_CATEGORY_B = 'B';
+const char ICAO_AIRCRAFT_CATEGORY_C = 'C';
+const char ICAO_AIRCRAFT_CATEGORY_D = 'D';
+const char ICAO_AIRCRAFT_CATEGORY_E = 'E';
+
+/**
+ * Calculate flight parameter based on aircraft performance data. 
+ * This is based on simple rules: it does not (yet) include data
+ * such as winds aloft, payload or temperature impact on engine
+ * performance. */
+class AircraftPerformance
+{
+public:
+    AircraftPerformance();
+
+    double turnRateDegSec() const;
+
+    double turnRadiusMForAltitude(int altitudeFt) const;
+    
+    double groundSpeedForAltitudeKnots(int altitudeFt) const;
+
+    int computePreviousAltitude(double distanceM, int targetAltFt) const;
+    int computeNextAltitude(double distanceM, int initialAltFt) const;
+
+    double distanceNmBetween(int initialElevationFt, int targetElevationFt) const;
+    double timeBetween(int initialElevationFt, int targetElevationFt) const;
+
+    double timeToCruise(double cruiseDistanceNm, int cruiseAltitudeFt) const;
+
+private:
+    void readPerformanceData();
+
+    void icaoCategoryData();
+
+    /**
+     * @brief heuristicCatergoryFromTags - based on the aircraft tags, figure
+     * out a plausible ICAO category. Returns cat A if nothing better could
+     * be determined.
+     * @return a string containing a single ICAO category character A..E
+     */
+    std::string heuristicCatergoryFromTags() const;
+
+    class Bracket
+    {
+    public:
+        Bracket(int atOrBelow, int climb, int descent, double speed, bool isMach = false) :
+            atOrBelowAltitudeFt(atOrBelow),
+            climbRateFPM(climb),
+            descentRateFPM(descent),
+            speedIASOrMach(speed),
+            speedIsMach(isMach)
+        { }
+        
+        int gsForAltitude(int altitude) const;
+
+        double climbTime(int alt1, int alt2) const;
+        double climbDistanceM(int alt1, int alt2) const;
+        double descendTime(int alt1, int alt2) const;
+        double descendDistanceM(int alt1, int alt2) const;
+
+        int atOrBelowAltitudeFt;
+        int climbRateFPM;
+        int descentRateFPM;
+        double speedIASOrMach;
+        bool speedIsMach = false;
+    };
+
+    using PerformanceVec = std::vector<Bracket>;
+
+    using BracketRange = std::pair<PerformanceVec::const_iterator, PerformanceVec::const_iterator>;
+
+    PerformanceVec::const_iterator bracketForAltitude(int altitude) const;
+    BracketRange rangeForAltitude(int lowAltitude, int highAltitude) const;
+
+
+    using TraversalFunc = std::function<void(const Bracket& bk, int alt1, int alt2)>;
+    void traverseAltitudeRange(int initialElevationFt, int targetElevationFt, TraversalFunc tf) const;
+
+
+    PerformanceVec _perfData;
+};
+
+}
+
+#endif // AIRCRAFTPERFORMANCE_HXX

src/Aircraft/CMakeLists.txt

@@ -6,6 +6,7 @@ set(SOURCES
 	flightrecorder.cxx
     FlightHistory.cxx
 		initialstate.cxx
+	AircraftPerformance.cxx
 	)
 
 set(HEADERS
@@ -14,6 +15,7 @@ set(HEADERS
 	flightrecorder.hxx
     FlightHistory.hxx
 		initialstate.hxx
+	AircraftPerformance.hxx
 	)
 
 

src/Navaids/FlightPlan.cxx

@@ -49,6 +49,7 @@
 #include <Navaids/routePath.hxx>
 #include <Navaids/airways.hxx>
 #include <Autopilot/route_mgr.hxx>
+#include <Aircraft/AircraftPerformance.hxx>
 
 using std::string;
 using std::vector;

src/Navaids/FlightPlan.hxx

@@ -35,12 +35,6 @@ class FlightPlan;
 
 typedef SGSharedPtr<FlightPlan> FlightPlanRef;
 
-const char ICAO_AIRCRAFT_CATEGORY_A = 'A';
-const char ICAO_AIRCRAFT_CATEGORY_B = 'B';
-const char ICAO_AIRCRAFT_CATEGORY_C = 'C';
-const char ICAO_AIRCRAFT_CATEGORY_D = 'D';
-const char ICAO_AIRCRAFT_CATEGORY_E = 'E';
-
 class FlightPlan : public RouteBase
 {
 public:

src/Navaids/routePath.cxx

@@ -1,7 +1,21 @@
+// routePath.hxx - compute data about planned route
+//
+// Copyright (C) 2018  James Turner  <james@flightgear.org>
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of the
+// License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 
-#ifdef HAVE_CONFIG_H
+#include "config.h"
-#  include "config.h"
-#endif
 
 #include <algorithm>
 
@@ -16,6 +30,7 @@
 #include <Navaids/waypoint.hxx>
 #include <Navaids/FlightPlan.hxx>
 #include <Navaids/positioned.hxx>
+#include <Aircraft/AircraftPerformance.hxx>
 
 namespace flightgear {
 
@@ -589,26 +604,6 @@ public:
 
 typedef std::vector<WayptData> WayptDataVec;
 
-class PerformanceBracket
-{
-public:
-  PerformanceBracket(double atOrBelow, double climb, double descent, double speed, bool isMach = false) :
-    atOrBelowAltitudeFt(atOrBelow),
-    climbRateFPM(climb),
-    descentRateFPM(descent),
-    speedIASOrMach(speed),
-    speedIsMach(isMach)
-  { }
-  
-  double atOrBelowAltitudeFt;
-  double climbRateFPM;
-  double descentRateFPM;
-  double speedIASOrMach;
-  bool speedIsMach;
-};
-
-typedef std::vector<PerformanceBracket> PerformanceBracketVec;
-
 bool isDescentWaypoint(const WayptRef& wpt)
 {
   return (wpt->flag(WPT_APPROACH) && !wpt->flag(WPT_MISS)) || wpt->flag(WPT_ARRIVAL);
@@ -619,28 +614,9 @@ class RoutePath::RoutePathPrivate
 public:
     WayptDataVec waypoints;
 
-    char aircraftCategory;
+    AircraftPerformance perf;
-    PerformanceBracketVec perf;
-    double pathTurnRate;
     bool constrainLegCourses;
 
-  PerformanceBracketVec::const_iterator
-  findPerformanceBracket(double altFt) const
-  {
-    PerformanceBracketVec::const_iterator r;
-    PerformanceBracketVec::const_iterator result = perf.begin();
-
-    for (r = perf.begin(); r != perf.end(); ++r) {
-      if (r->atOrBelowAltitudeFt > altFt) {
-        break;
-      }
-      
-      result = r;
-    } // of brackets iteration
-    
-    return result;
-  }
-  
   void computeDynamicPosition(int index)
   {
     auto previous(previousValidWaypoint(index));
@@ -653,18 +629,7 @@ public:
       HeadingToAltitude* h = (HeadingToAltitude*) wpt.get();
       
       double altFt = computeVNAVAltitudeFt(index - 1);
-      double altChange = h->altitudeFt() - altFt;
+      double distanceM = perf.distanceNmBetween(altFt, h->altitudeFt()) * SG_NM_TO_METER;
-      PerformanceBracketVec::const_iterator it = findPerformanceBracket(altFt);
-      double speedMSec = groundSpeedForAltitude(altFt) * SG_KT_TO_MPS;
-      double timeToChangeSec;
-      
-      if (isDescentWaypoint(wpt)) {
-        timeToChangeSec = (altChange / it->descentRateFPM) * 60.0;
-      } else {
-        timeToChangeSec = (altChange / it->climbRateFPM) * 60.0;
-      }
-
-      double distanceM = timeToChangeSec * speedMSec;
       double hdg = h->headingDegMagnetic() + magVarFor(previous->pos);
       waypoints[index].pos = SGGeodesy::direct(previous->turnExitPos, hdg, distanceM);
       waypoints[index].posValid = true;
@@ -742,12 +707,7 @@ public:
       
       double fixedAlt = altitudeForIndex(next);
       double distanceM = distanceBetweenIndices(index, next);
-      double speedMSec = groundSpeedForAltitude(fixedAlt) * SG_KT_TO_MPS;
+      return perf.computePreviousAltitude(distanceM, fixedAlt);
-      double minutes = (distanceM / speedMSec) / 60.0;
-      
-      PerformanceBracketVec::const_iterator it = findPerformanceBracket(fixedAlt);
-      return fixedAlt + (it->descentRateFPM * minutes);
-
     } else {
       // climb
       int prev = findPreceedingKnownAltitude(index);
@@ -757,13 +717,8 @@ public:
       
       double fixedAlt = altitudeForIndex(prev);
       double distanceM = distanceBetweenIndices(prev, index);
-      double speedMSec = groundSpeedForAltitude(fixedAlt) * SG_KT_TO_MPS;
+      return perf.computeNextAltitude(distanceM, fixedAlt);
-      double minutes = (distanceM / speedMSec) / 60.0;
-      
-      PerformanceBracketVec::const_iterator it = findPerformanceBracket(fixedAlt);
-      return fixedAlt + (it->climbRateFPM * minutes);
     }
-    
   }
   
   int findPreceedingKnownAltitude(int index) const
@@ -831,47 +786,6 @@ public:
     return 0.0;
   }
   
-  double groundSpeedForAltitude(double altitude) const
-  {
-      PerformanceBracketVec::const_iterator it = findPerformanceBracket(altitude);
-      if (it->speedIsMach) {
-          return 300.0;
-      } else {
-          // FIXME - convert IAS to ground-speed, using standard atmosphere / temperature model
-          return it->speedIASOrMach;
-      }
-
-#if 0
-    if (0) {
-      double mach;
-      
-      if (it->speedIsMach) {
-        mach = it->speedIASOrMach; // easy
-      } else {
-        const double Cs_0 = 661.4786; // speed of sound at sea level, knots
-        const double P_0 = 29.92126;
-        const double P = P_0 * pow(, );
-        // convert IAS (which we will treat as CAS) to Mach based on altitude
-      }
-      
-      double oatK;
-      double Cs = sqrt(SG_gamma * SG_R_m2_p_s2_p_K * oatK);
-      
-      double tas = mach * Cs;
-      
-/*
-      P_0= 29.92126 "Hg = 1013.25 mB = 2116.2166 lbs/ft^2
-      P= P_0*(1-6.8755856*10^-6*PA)^5.2558797, pressure altitude, PA<36,089.24ft
-      CS= 38.967854*sqrt(T+273.15)  where T is the (static/true) OAT in Celsius.
-      
-      DP=P_0*((1 + 0.2*(IAS/CS_0)^2)^3.5 -1)
-      M=(5*( (DP/P + 1)^(2/7) -1) )^0.5   (*)
-      TAS= M*CS
-*/
-    }
-#endif
-  }
-
   double distanceBetweenIndices(int from, int to) const
   {
     double total = 0.0;
@@ -883,40 +797,6 @@ public:
     return total;
   }
   
-  void initPerfData()
-  {
-      pathTurnRate = 3.0; // 3 deg/sec = 180deg/min = standard rate turn
-      switch (aircraftCategory) {
-      case ICAO_AIRCRAFT_CATEGORY_A:
-          perf.push_back(PerformanceBracket(4000, 600, 1200, 75));
-          perf.push_back(PerformanceBracket(10000, 600, 1200, 140));
-          break;
-
-      case ICAO_AIRCRAFT_CATEGORY_B:
-          perf.push_back(PerformanceBracket(4000, 100, 1200, 100));
-          perf.push_back(PerformanceBracket(10000, 800, 1200, 160));
-          perf.push_back(PerformanceBracket(18000, 600, 1800, 200));
-          break;
-
-      case ICAO_AIRCRAFT_CATEGORY_C:
-          perf.push_back(PerformanceBracket(4000, 1800, 1800, 150));
-          perf.push_back(PerformanceBracket(10000, 1800, 1800, 200));
-          perf.push_back(PerformanceBracket(18000, 1200, 1800, 270));
-          perf.push_back(PerformanceBracket(60000, 800, 1200, 0.80, true /* is Mach */));
-          break;
-
-      case ICAO_AIRCRAFT_CATEGORY_D:
-      case ICAO_AIRCRAFT_CATEGORY_E:
-      default:
-
-          perf.push_back(PerformanceBracket(4000, 1800, 1800, 180));
-          perf.push_back(PerformanceBracket(10000, 1800, 1800, 230));
-          perf.push_back(PerformanceBracket(18000, 1200, 1800, 270));
-          perf.push_back(PerformanceBracket(60000, 800, 1200, 0.87, true /* is Mach */));
-          break;
-      }
-  }
-
     WayptDataVec::iterator previousValidWaypoint(unsigned int index)
     {
         if (index == 0) {
@@ -968,7 +848,6 @@ RoutePath::RoutePath(const flightgear::FlightPlan* fp) :
         d->waypoints.push_back(WayptData(wpt));
     }
 
-    d->aircraftCategory = fp->icaoAircraftCategory()[0];
     d->constrainLegCourses = fp->followLegTrackToFixes();
     commonInit();
 }
@@ -979,15 +858,13 @@ RoutePath::~RoutePath()
 
 void RoutePath::commonInit()
 {  
-  d->initPerfData();
-  
   WayptDataVec::iterator it;
   for (it = d->waypoints.begin(); it != d->waypoints.end(); ++it) {
     it->initPass0();
   }
   
   for (unsigned int i=1; i<d->waypoints.size(); ++i) {
-    WayptData* nextPtr = ((i + 1) < d->waypoints.size()) ? &d->waypoints[i+1] : 0;
+    WayptData* nextPtr = ((i + 1) < d->waypoints.size()) ? &d->waypoints[i+1] : nullptr;
     d->waypoints[i].initPass1(d->waypoints[i-1], nextPtr);
   }
 
@@ -997,8 +874,7 @@ void RoutePath::commonInit()
       }
 
       double alt = 0.0; // FIXME
-      double gs = d->groundSpeedForAltitude(alt);
+      double radiusM = d->perf.turnRadiusMForAltitude(alt);
-      double radiusM = ((360.0 / d->pathTurnRate) * gs * SG_KT_TO_MPS) / SGMiscd::twopi();
 
       if (i > 0) {
           auto prevIt = d->previousValidWaypoint(i);
@@ -1142,12 +1018,12 @@ SGGeodVec RoutePath::pathForHold(Hold* hold) const
     double hdg = hold->inboundRadial();
     double turnDelta = 180.0 / turnSteps;
     double altFt = 0.0; // FIXME
-    double gsKts = d->groundSpeedForAltitude(altFt);
+    double gsKts = d->perf.groundSpeedForAltitudeKnots(altFt);
 
     SGGeodVec r;
     double az2;
-    double stepTime = turnDelta / d->pathTurnRate; // in seconds
+    double stepTime = turnDelta / 3.0; // 3.0 is degrees/sec for standard rate turn
-    double stepDist = gsKts * (stepTime / 3600.0) * SG_NM_TO_METER;
+    double stepDist = gsKts * SG_KT_TO_MPS * stepTime;
     double legDist = hold->isDistance() ?
     hold->timeOrDistance()
     : gsKts * (hold->timeOrDistance() / 3600.0);

test_suite/unit_tests/Navaids/CMakeLists.txt

@@ -3,11 +3,13 @@ set(TESTSUITE_SOURCES
     ${CMAKE_CURRENT_SOURCE_DIR}/TestSuite.cxx
     ${CMAKE_CURRENT_SOURCE_DIR}/test_flightplan.cxx
     ${CMAKE_CURRENT_SOURCE_DIR}/test_navaids2.cxx
+    ${CMAKE_CURRENT_SOURCE_DIR}/test_aircraftPerformance.cxx
     PARENT_SCOPE
 )
 
 set(TESTSUITE_HEADERS
     ${TESTSUITE_HEADERS}
     ${CMAKE_CURRENT_SOURCE_DIR}/test_flightplan.hxx
+    ${CMAKE_CURRENT_SOURCE_DIR}/test_aircraftPerformance.hxx
     PARENT_SCOPE
 )

test_suite/unit_tests/Navaids/TestSuite.cxx

@@ -19,8 +19,9 @@
 
 #include "test_flightplan.hxx"
 #include "test_navaids2.hxx"
-
+#include "test_aircraftPerformance.hxx"
 
 // Set up the unit tests.
 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION(FlightplanTests, "Unit tests");
 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION(NavaidsTests, "Unit tests");
+CPPUNIT_TEST_SUITE_NAMED_REGISTRATION(AircraftPerformanceTests, "Unit tests");

test_suite/unit_tests/Navaids/test_aircraftPerformance.cxx

@@ -0,0 +1,73 @@
+#include "test_aircraftPerformance.hxx"
+
+#include "test_suite/FGTestApi/globals.hxx"
+#include "test_suite/FGTestApi/NavDataCache.hxx"
+
+#include <simgear/misc/strutils.hxx>
+
+#include <Main/fg_props.hxx>
+#include <Aircraft/AircraftPerformance.hxx>
+
+using namespace flightgear;
+
+// Set up function for each test.
+void AircraftPerformanceTests::setUp()
+{
+    FGTestApi::setUp::initTestGlobals("aircraft-perf");
+    FGTestApi::setUp::initNavDataCache();
+}
+
+
+// Clean up after each test.
+void AircraftPerformanceTests::tearDown()
+{
+    FGTestApi::tearDown::shutdownTestGlobals();
+}
+
+void AircraftPerformanceTests::testBasic()
+{
+    fgSetString("/aircraft/performance/icao-category", "C");
+    AircraftPerformance ap;
+
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.groundSpeedForAltitudeKnots(1000), 152, 1e-3);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.groundSpeedForAltitudeKnots(8000), 224, 1e-3);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.groundSpeedForAltitudeKnots(20000), 491, 1e-3);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.groundSpeedForAltitudeKnots(35000), 461, 1e-3);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.groundSpeedForAltitudeKnots(38000), 458, 1e-3);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.groundSpeedForAltitudeKnots(40000), 458, 1e-3);
+
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.timeBetween(3000, 6000), 100, 1);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.distanceNmBetween(3000, 6000), 5.430, 1e-3);
+    
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.timeBetween(36000, 34000), 100, 1);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.distanceNmBetween(36000, 34000), 12.805, 1e-1);
+    
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.timeBetween(15000, 20000), 300, 1);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.distanceNmBetween(15000, 18000), 14.270, 1);
+
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.timeBetween(2000, 25000), 1191.6, 1);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.distanceNmBetween(2000, 25000), 123.06, 1);
+    
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.timeBetween(36000, 3000), 1666.6, 1);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(ap.distanceNmBetween(36000, 3000), 162.02, 1);
+    
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(251.5, ap.timeToCruise(32.0, 350000), 1);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(503.0, ap.timeToCruise(64.0, 380000), 1);
+
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(1553.7, ap.turnRadiusMForAltitude(4000), 1);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(3322.4, ap.turnRadiusMForAltitude(16000), 1);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(4602.6, ap.turnRadiusMForAltitude(30000), 1);
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(4475.5, ap.turnRadiusMForAltitude(38000), 1);
+}
+
+void AircraftPerformanceTests::testAltitudeGradient()
+{
+    fgSetString("/aircraft/performance/icao-category", "E");
+    AircraftPerformance ap;
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(8332, ap.computePreviousAltitude(10000, 6000), 1);
+    
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(3260, ap.computeNextAltitude(4000, 2000), 1);
+
+
+}
+

test_suite/unit_tests/Navaids/test_aircraftPerformance.hxx

@@ -0,0 +1,49 @@
+/*
+ * Copyright (C) 2018 James Turner
+ *
+ * This file is part of the program FlightGear.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+
+#ifndef _FG_AIRCRAFT_PERFORMANCE_UNIT_TESTS_HXX
+#define _FG_AIRCRAFT_PERFORMANCE_UNIT_TESTS_HXX
+
+
+#include <cppunit/extensions/HelperMacros.h>
+#include <cppunit/TestFixture.h>
+
+
+class AircraftPerformanceTests : public CppUnit::TestFixture
+{
+    // Set up the test suite.
+    CPPUNIT_TEST_SUITE(AircraftPerformanceTests);
+    CPPUNIT_TEST(testBasic);
+    CPPUNIT_TEST(testAltitudeGradient);
+    CPPUNIT_TEST_SUITE_END();
+
+public:
+    // Set up function for each test.
+    void setUp();
+
+    // Clean up after each test.
+    void tearDown();
+
+    // The tests.
+    void testBasic();
+    void testAltitudeGradient();
+};
+
+#endif  // AircraftPerformanceTests