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.Altitudes = None
        self.Directions = None
        self.Windspeeds = None
        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.Weather or None == self.Weather.Provider:
            return

        if None == self.LastWeatherUpdate or self.LastWeatherUpdate != self.Weather.Provider.UpdateTime:
            self.MinimumAltitude = 1000000
            self.MaximumAltitude = -1
            self.WindDirectionModel = None
            self.WindSpeedModel = None
            self.Altitudes = None
            self.Directions = None
            self.Windspeeds = None

            if None != self.Weather.Provider.WindData and self.Gafor in self.Weather.Provider.WindData:
                self.Altitudes = []
                self.Directions = []
                self.Windspeeds = []

                # collect the data for the wind model
                for level in self.Weather.Provider.WindData[self.Gafor]:
                    self.Altitudes.append(level[0])
                    self.Directions.append(level[1])
                    self.Windspeeds.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(self.Altitudes):
                    self.WindDirectionModel = scipy.interpolate.interp1d(self.Altitudes, self.Directions)
                    self.WindSpeedModel = scipy.interpolate.interp1d(self.Altitudes, self.Windspeeds)
                    self.LastWeatherUpdate = self.Weather.Provider.UpdateTime
            else:
                self.LastWeatherUpdate = None

    def interpolateWindData(self, altitude : int):
        self.updateWindModel()

        # initialized 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

        return speed, direction

    def calculateGS(self, altitude : int, ias : int, heading : int):
        speed, direction = self.interpolateWindData(altitude)
        tas = self.calculateTAS(altitude, ias)
        return tas + speed * math.cos(math.radians(direction) - math.radians(heading))

    def convertGSToTAS(self, altitude : int, gs : int, heading : int):
        speed, direction = self.interpolateWindData(altitude)
        return gs - speed * math.cos(math.radians(direction) - math.radians(heading))

    def estimateCourse(self, altitude : int, gs : int, heading : int):
        tas = self.convertGSToTAS(altitude, gs, heading)
        speed, direction = self.interpolateWindData(altitude)

        aca = heading - direction
        wca = speed * aca / tas

        if 0 <= aca:
            course = heading + wca
        else:
            course = heading - wca

        while 0 > course:
            course += 360
        while 360 < course:
            course -= 360

        return course