aman-sys/aman/sys/WeatherModel.py

#!/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))