From e3d032942e9edb87c15e6c79c56cfcdf6364e31f Mon Sep 17 00:00:00 2001 From: James Turner Date: Thu, 16 Aug 2018 19:43:09 -0700 Subject: [PATCH] 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 | 390 ++++++++++++++++++ src/Aircraft/AircraftPerformance.hxx | 112 +++++ src/Aircraft/CMakeLists.txt | 2 + src/Navaids/FlightPlan.cxx | 1 + src/Navaids/FlightPlan.hxx | 6 - src/Navaids/routePath.cxx | 180 ++------ test_suite/unit_tests/Navaids/CMakeLists.txt | 2 + test_suite/unit_tests/Navaids/TestSuite.cxx | 3 +- .../Navaids/test_aircraftPerformance.cxx | 73 ++++ .../Navaids/test_aircraftPerformance.hxx | 49 +++ 10 files changed, 659 insertions(+), 159 deletions(-) create mode 100644 src/Aircraft/AircraftPerformance.cxx create mode 100644 src/Aircraft/AircraftPerformance.hxx create mode 100644 test_suite/unit_tests/Navaids/test_aircraftPerformance.cxx create mode 100644 test_suite/unit_tests/Navaids/test_aircraftPerformance.hxx diff --git a/src/Aircraft/AircraftPerformance.cxx b/src/Aircraft/AircraftPerformance.cxx new file mode 100644 index 000000000..364382307 --- /dev/null +++ b/src/Aircraft/AircraftPerformance.cxx @@ -0,0 +1,390 @@ +// AircraftPerformance.cxx - compute data about planned acft performance +// +// Copyright (C) 2018 James Turner +// 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 + +#include + +#include
+ +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(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(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; +} + diff --git a/src/Aircraft/AircraftPerformance.hxx b/src/Aircraft/AircraftPerformance.hxx new file mode 100644 index 000000000..c2183de8a --- /dev/null +++ b/src/Aircraft/AircraftPerformance.hxx @@ -0,0 +1,112 @@ +// AircraftPerformance.hxx - compute data about planned acft performance +// +// Copyright (C) 2018 James Turner +// 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 +#include + +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; + + using BracketRange = std::pair; + + PerformanceVec::const_iterator bracketForAltitude(int altitude) const; + BracketRange rangeForAltitude(int lowAltitude, int highAltitude) const; + + + using TraversalFunc = std::function; + void traverseAltitudeRange(int initialElevationFt, int targetElevationFt, TraversalFunc tf) const; + + + PerformanceVec _perfData; +}; + +} + +#endif // AIRCRAFTPERFORMANCE_HXX diff --git a/src/Aircraft/CMakeLists.txt b/src/Aircraft/CMakeLists.txt index 419923663..2038dd4d7 100644 --- a/src/Aircraft/CMakeLists.txt +++ b/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 ) diff --git a/src/Navaids/FlightPlan.cxx b/src/Navaids/FlightPlan.cxx index a745a9faf..6f32a73ee 100644 --- a/src/Navaids/FlightPlan.cxx +++ b/src/Navaids/FlightPlan.cxx @@ -49,6 +49,7 @@ #include #include #include +#include using std::string; using std::vector; diff --git a/src/Navaids/FlightPlan.hxx b/src/Navaids/FlightPlan.hxx index b8d58e596..4bee2941b 100644 --- a/src/Navaids/FlightPlan.hxx +++ b/src/Navaids/FlightPlan.hxx @@ -35,12 +35,6 @@ class FlightPlan; typedef SGSharedPtr 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: diff --git a/src/Navaids/routePath.cxx b/src/Navaids/routePath.cxx index 87779c429..a0749da92 100644 --- a/src/Navaids/routePath.cxx +++ b/src/Navaids/routePath.cxx @@ -1,7 +1,21 @@ +// routePath.hxx - compute data about planned route +// +// Copyright (C) 2018 James Turner +// 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" -#endif +#include "config.h" #include @@ -16,6 +30,7 @@ #include #include #include +#include namespace flightgear { @@ -589,26 +604,6 @@ public: typedef std::vector 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 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; - PerformanceBracketVec perf; - double pathTurnRate; + AircraftPerformance perf; 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; - 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 distanceM = perf.distanceNmBetween(altFt, h->altitudeFt()) * SG_NM_TO_METER; 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; - double minutes = (distanceM / speedMSec) / 60.0; - - PerformanceBracketVec::const_iterator it = findPerformanceBracket(fixedAlt); - return fixedAlt + (it->descentRateFPM * minutes); - + return perf.computePreviousAltitude(distanceM, fixedAlt); } 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; - double minutes = (distanceM / speedMSec) / 60.0; - - PerformanceBracketVec::const_iterator it = findPerformanceBracket(fixedAlt); - return fixedAlt + (it->climbRateFPM * minutes); + return perf.computeNextAltitude(distanceM, fixedAlt); } - } 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(); } @@ -978,16 +857,14 @@ 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; iwaypoints.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 = ((360.0 / d->pathTurnRate) * gs * SG_KT_TO_MPS) / SGMiscd::twopi(); + double radiusM = d->perf.turnRadiusMForAltitude(alt); 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 stepDist = gsKts * (stepTime / 3600.0) * SG_NM_TO_METER; + double stepTime = turnDelta / 3.0; // 3.0 is degrees/sec for standard rate turn + double stepDist = gsKts * SG_KT_TO_MPS * stepTime; double legDist = hold->isDistance() ? hold->timeOrDistance() : gsKts * (hold->timeOrDistance() / 3600.0); diff --git a/test_suite/unit_tests/Navaids/CMakeLists.txt b/test_suite/unit_tests/Navaids/CMakeLists.txt index fde0718eb..bf5556a73 100644 --- a/test_suite/unit_tests/Navaids/CMakeLists.txt +++ b/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 ) diff --git a/test_suite/unit_tests/Navaids/TestSuite.cxx b/test_suite/unit_tests/Navaids/TestSuite.cxx index b024ff6c7..e4a1f93a5 100644 --- a/test_suite/unit_tests/Navaids/TestSuite.cxx +++ b/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"); diff --git a/test_suite/unit_tests/Navaids/test_aircraftPerformance.cxx b/test_suite/unit_tests/Navaids/test_aircraftPerformance.cxx new file mode 100644 index 000000000..871eeea2d --- /dev/null +++ b/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 + +#include
+#include + +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); + + +} + diff --git a/test_suite/unit_tests/Navaids/test_aircraftPerformance.hxx b/test_suite/unit_tests/Navaids/test_aircraftPerformance.hxx new file mode 100644 index 000000000..541a7460b --- /dev/null +++ b/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 . + */ + + +#ifndef _FG_AIRCRAFT_PERFORMANCE_UNIT_TESTS_HXX +#define _FG_AIRCRAFT_PERFORMANCE_UNIT_TESTS_HXX + + +#include +#include + + +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