From a86dfa01d867156a66d43ddbe08472f2a23256cd Mon Sep 17 00:00:00 2001
From: Sven Czarnian <devel@svcz.de>
Date: Mon, 11 Oct 2021 21:58:14 +0200
Subject: [PATCH] add a weather model

---
 aman/sys/WeatherModel.py | 148 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 148 insertions(+)
 create mode 100644 aman/sys/WeatherModel.py

diff --git a/aman/sys/WeatherModel.py b/aman/sys/WeatherModel.py
new file mode 100644
index 0000000..32d4904
--- /dev/null
+++ b/aman/sys/WeatherModel.py
@@ -0,0 +1,148 @@
+#!/usr/bin/env python
+
+from aman.com.Weather import Weather
+
+import math
+import scipy.interpolate
+
+class WeatherModel:
+    def __init__(self, gaforId, weather : Weather):
+        self.gafor = gaforId
+        self.weather = weather
+        self.windDirectionModel = None
+        self.windSpeedModel = None
+        self.lastWeatherUpdate = None
+        self.minimumAltitude = 1000000
+        self.maximumAltitude = -1
+
+        # create the density interpolation model
+        # the density model is based on https://aerotoolbox.com/atmcalc/
+        altitudes = [
+            50000,
+            45000,
+            40000,
+            38000,
+            36000,
+            34000,
+            32000,
+            30000,
+            28000,
+            26000,
+            24000,
+            22000,
+            20000,
+            18000,
+            16000,
+            15000,
+            14000,
+            13000,
+            12000,
+            11000,
+            10000,
+            9000,
+            8000,
+            7000,
+            6000,
+            5000,
+            4000,
+            3000,
+            2000,
+            1000,
+            0
+        ]
+        densities = [
+            0.18648,
+            0.23714,
+            0.24617,
+            0.33199,
+            0.36518,
+            0.39444,
+            0.42546,
+            0.45831,
+            0.402506,
+            0.432497,
+            0.464169,
+            0.60954,
+            0.65269,
+            0.69815,
+            0.74598,
+            0.77082,
+            0.79628,
+            0.82238,
+            0.84914,
+            0.87655,
+            0.90464,
+            0.93341,
+            0.96287,
+            0.99304,
+            1.02393,
+            1.05555,
+            1.08791,
+            1.12102,
+            1.1549,
+            1.18955,
+            1.225
+        ]
+        self.densityModel = scipy.interpolate.interp1d(altitudes, densities)
+
+    def calculateTAS(self, altitude : int, ias : int):
+        if altitude >= 50000:
+            altitude = 49999
+        if altitude <= 0:
+            altitude = 1
+
+        # calculation based on https://aerotoolbox.com/airspeed-conversions/
+        return ias * math.sqrt(1.225 / self.densityModel(altitude).item())
+
+    def updateWindModel(self):
+        if None == self.lastWeatherUpdate or self.lastWeatherUpdate != self.weather.provider.updateTime:
+            self.lastWeatherUpdate = self.weather.provider.updateTime
+
+            self.minimumAltitude = 1000000
+            self.maximumAltitude = -1
+            self.windDirectionModel = None
+            self.windSpeedModel = None
+
+            if None != self.weather.provider.windData and self.gafor in self.weather.provider.windData:
+                altitudes = []
+                directions = []
+                speeds = []
+
+                # collect the data for the wind model
+                for level in self.weather.provider.windData[self.gafor]:
+                    altitudes.append(level[0])
+                    directions.append(level[1])
+                    speeds.append(level[2])
+
+                    # define the thresholds for later boundary checks
+                    if self.minimumAltitude > level[0]:
+                        self.minimumAltitude = level[0]
+                    if self.maximumAltitude < level[0]:
+                        self.maximumAltitude = level[0]
+
+                # calculate the models
+                if 1 < len(altitudes):
+                    self.windDirectionModel = scipy.interpolate.interp1d(altitudes, directions)
+                    self.windSpeedModel = scipy.interpolate.interp1d(altitudes, speeds)
+
+    def calculateGS(self, altitude : int, ias : int, heading : int):
+        self.updateWindModel()
+        tas = self.calculateTAS(altitude, ias)
+
+        # initialize the wind data
+        if None != self.windDirectionModel and None != self.windSpeedModel:
+            direction = 0.0
+            speed = 0.0
+            if None != self.windSpeedModel and None != self.windDirectionModel:
+                if self.maximumAltitude <= altitude:
+                    altitude = self.maximumAltitude - 1
+                if self.minimumAltitude >= altitude:
+                    altitude = self.minimumAltitude + 1
+                direction = self.windDirectionModel(altitude).item()
+                speed = self.windSpeedModel(altitude).item()
+        else:
+            speed = 0
+            direction = 0
+
+        # calculate the ground speed based on the headwind component
+        return tas + speed * math.cos(math.radians(direction) - math.radians(heading))