Speed of sound, viscosity and custom atmosphere update

README.md

@@ -1,9 +1,10 @@
 # ISA Calculator
-Basic ISA Calculator for Python projects
+_Basic ISA Calculator for Python projects_
 
 Contents:
-  - Atmospheric Model up to 120km
+  - Atmospheric Model up to 110km
   - Accurate Calculations for Temperature, Pressure and Density
+  - Custom defined atmospheric conditions
 
 
 With this module, it is possible to calculate, using the 1976 standard atmosphere model, the Temperature,
@@ -16,21 +17,23 @@ To use this package, follow these steps:
     - Install isacalc
     - Import isacalc as isa
     - Define the Atmosphere Model: isa.get_atmosphere()
-    - Calculate the Temperature, Pressure and Density at height h using isa.calculate_at_h(h, atmosphere_model)
+    - It calculates, at the defined height:
+        - Temperature
+        - Pressure
+        - Density
+        - Speed of sound
+        - Dynamic Viscosity
 
 And thats it! An example Script:
 
 
     import isacalc as isa
 
     atmosphere = isa.get_atmosphere()
+
     h = 50000.0
 
-    h, T, P, d = isa.calculate_at_h(h, atmosphere)
-
+    T, P, d, a, mu = isa.calculate_at_h(h, atmosphere)
 
-Planned Future Features:
 
-    - Atmosphere Model Customization
-    - Mach number and viscosity calculations
 

isacalc/main_executable.py

@@ -1,11 +1,12 @@
 from isacalc.src import Atmosphere
 
-def get_atmosphere() -> Atmosphere:
+def get_atmosphere(*args, **kwargs) -> Atmosphere:
     """
     Function to obtain the atmosphere model
+    :param kwargs: To define a custom atmosphere model
     :return: Atmospheric Model
     """
-    return Atmosphere()
+    return Atmosphere(*args, **kwargs)
 
 
 def calculate_at_h(h: float, atmosphere_model: Atmosphere = get_atmosphere()) -> list:
@@ -21,6 +22,6 @@ def calculate_at_h(h: float, atmosphere_model: Atmosphere = get_atmosphere()) ->
 if __name__ == "__main__":
 
     atmosphere = get_atmosphere()
-    h = 50000.0
+    h = 80000
 
-    h, T, P, d = calculate_at_h(h, atmosphere)
+    T, P, d, a, mu = calculate_at_h(h, atmosphere)

isacalc/src/atmosphere.py

@@ -4,19 +4,53 @@
 
 class Atmosphere(object):
 
-    def __init__(self):
+    def __init__(self, *args, **kwargs):
 
-        self.__p0 = 101325.0
-        self.__d0 = 1.225
+        if kwargs:
+            # User Defined Atmosphere Model
+            self.__p0 = kwargs['p0']
+            self.__d0 = kwargs['d0']
 
-        self.__Nn = np.array(["Troposphere", "Tropopause", "Stratosphere", "Stratosphere", "Stratopause", "Mesosphere", "Mesosphere", "Mesopause", "Thermosphere", "Thermosphere", "Thermosphere"])
-        self.__Tn = np.array([288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65, 186.95, 186.95, 201.95, 251.95])
-        self.__Hn = np.array([0, 11000.0, 20000.0, 32000.0, 47000.0, 51000.0, 71000.0, 84852.0, 90000.0, 100000.0, 110000.0])
-        self.__Lt = np.array([1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1])
+            self.__Hn = kwargs['heights']
+            self.__Tn = kwargs['temp']
+
+            try:
+                self.__Nn = kwargs['names']
+
+            except KeyError:
+                self.__Nn = ['Noname']*len(self.__Hn)
+
+        else:
+            # Standard Values
+            self.__p0 = 101325.0
+            self.__d0 = 1.225
+
+            self.__Nn = ["Troposphere", "Tropopause", "Stratosphere", "Stratosphere", "Stratopause", "Mesosphere", "Mesosphere", "Mesopause", "Thermosphere", "Thermosphere", "Thermosphere"]
+            self.__Tn = [288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65, 186.95, 186.95, 201.95, 251.95]
+            self.__Hn = [0, 11000.0, 20000.0, 32000.0, 47000.0, 51000.0, 71000.0, 84852.0, 90000.0, 100000.0, 110000.0]
+
+        self.__Lt = self.__get_lapse(self.__Hn, self.__Tn)
 
         self.__layers = []
         self.__build()
 
+
+    @staticmethod
+    def __get_lapse(Hn, Tn) -> np.array:
+
+        types = []
+
+        for i in range(len(Hn)-1):
+            lapse = (Tn[i+1] - Tn[i])/(Hn[i+1] - Hn[i])
+
+            if lapse != 0:
+                types.append(1)
+
+            elif lapse == 0:
+                types.append(0)
+
+        return types
+
     def __build(self) -> None:
 
         p0, d0 = self.__p0, self.__d0
@@ -33,7 +67,7 @@ def __build(self) -> None:
                                     top_temperature=T_i,
                                     name=name)
 
-                _, _, p0, d0 = Layer.get_ceiling_values()
+                T0, p0, d0, a0, mu0 = Layer.get_ceiling_values()
 
             elif layer_type == 0:
 
@@ -44,7 +78,7 @@ def __build(self) -> None:
                                         max_height=h_i,
                                         name=name)
 
-                _, _, p0, d0 = Layer.get_ceiling_values()
+                T0, p0, d0, a0, mu0 = Layer.get_ceiling_values()
 
             else:
                 raise ValueError
@@ -53,9 +87,15 @@ def __build(self) -> None:
 
     def calculate(self, h) -> list:
 
+        if h > self.__Hn[-1] or h < self.__Hn[0]:
+            raise ValueError("Height is out of bounds")
+
         for idx in range(len(self.__layers)):
 
-            if self.__Hn[idx] <= h < self.__Hn[idx + 1]:
+            if h == self.__Hn[idx+1]:
+                return self.__layers[idx].get_ceiling_values()
+
+            elif self.__Hn[idx] < h < self.__Hn[idx + 1]:
                 return self.__layers[idx].get_intermediate_values(h)
 
             elif h >= self.__Hn[idx + 1]:

isacalc/src/layers.py

@@ -20,6 +20,7 @@ def __init__(self, base_height: float, base_temperature: float, base_pressure: f
 
         self.g0 = 9.80665
         self.R = 287.0
+        self.gamma = 1.4
 
         self.__h0 = base_height
         self.__T0 = base_temperature
@@ -29,12 +30,33 @@ def __init__(self, base_height: float, base_temperature: float, base_pressure: f
 
         self.__name = name
 
+    @staticmethod
+    def sutherland_viscosity(T, mu0=1.716e-5, T0=273.15, S=110.4):
+        """
+        Method to calculate the dynamic viscosity of air
+        :param T:   Temperature at which to calculate
+        :param mu0: Reference viscosity
+        :param T0:  Reference temperature
+        :param S:   Sutherland Temperature
+        :return:    Viscosity of air in kg/(m*s)
+        """
+        return mu0*(T/T0)**(1.5)*(T0+S)/(T+S)
+
+
+    def speed_of_sound(self, T):
+        """
+        Method to calculate the speed of sound at a certain temperature
+        :param T: Temperature in K
+        :return: speed of sound in m/s
+        """
+        return np.sqrt(self.gamma*self.R*T)
+
     def get_base_values(self) -> list:
         """
         Getter function to obtain the hidden layer states
         :return: List of all base values
         """
-        return [self.__h0, self.__T0, self.__p0, self.__d0]
+        return [self.__T0, self.__p0, self.__d0, self.speed_of_sound(self.__T0), self.sutherland_viscosity(self.__T0)]
 
     def get_ceiling_height(self) -> float:
         """
@@ -89,37 +111,45 @@ def __init__(self, base_height: float,
     def get_ceiling_values(self) -> list:
         """
         Method to get the temperature, pressure and density at the ceiling of the layer
-        :return:
+        :return: temperature, pressure, density, speed of sound
         """
 
-        h0, T0, P0, D0 = self.get_base_values()
+        T0, P0, D0, a0, mu0 = self.get_base_values()
         h = self.get_ceiling_height()
+        h0 = self.get_base_height()
 
         P = P0 * np.exp(-self.g0 / (self.R * T0) * (h - h0))
         D = D0 * np.exp(-self.g0 / (self.R * T0) * (h - h0))
 
-        return [h, T0, P, D]
+        a = self.speed_of_sound(T0)
+        mu = self.sutherland_viscosity(T0)
+
+        return [T0, P, D, a, mu]
 
     def get_intermediate_values(self, h) -> list:
         """
         Method to get the temperature, pressure and density at height h, between the base and ceiling of the layer
         :param h: Height at which to evaluate the temperature, pressure, density
-        :return:
+        :return: temperature, pressure, density, speed of sound
         """
         h_max = self.get_ceiling_height()
+        h0 = self.get_base_height()
 
         if h > h_max:
             raise ValueError
 
-        if h == self.get_base_height():
+        if h == h0:
             return self.get_base_values()
 
-        h0, T0, P0, D0 = self.get_base_values()
+        T0, P0, D0, a0, mu0 = self.get_base_values()
 
         P = P0 * np.exp(-self.g0 / (self.R * T0) * (h - h0))
         D = D0 * np.exp(-self.g0 / (self.R * T0) * (h - h0))
 
-        return [h, T0, P, D]
+        a = self.speed_of_sound(T0)
+        mu = self.sutherland_viscosity(T0)
+
+        return [T0, P, D, a, mu]
 
 
 class NormalLayer(Layer):
@@ -144,35 +174,41 @@ def __init__(self, base_height: float,
     def get_ceiling_values(self) -> list:
         """
         Method to get the temperature, pressure and density at the ceiling of the layer
-        :return:
+        :return: temperature, pressure, density, speed of sound
         """
 
-        h0, T0, P0, D0 = self.get_base_values()
+        T0, P0, D0, a0, mu0 = self.get_base_values()
         h = self.get_ceiling_height()
+        h0 = self.get_base_height()
 
         L = (self.__T_top - T0) / (h - h0)
         C = -self.g0 / (L * self.R)  # To Simplify and shorten code we define the following expression for the exponent
 
         P = P0 * (self.__T_top / T0) ** C
         D = D0 * (self.__T_top / T0) ** (C - 1)
 
-        return [h, self.__T_top, P, D]
+        a = self.speed_of_sound(self.__T_top)
+        mu = self.sutherland_viscosity(self.__T_top)
+
+        return [self.__T_top, P, D, a, mu]
 
     def get_intermediate_values(self, h) -> list:
         """
         Method to get the temperature, pressure and density at height h, between the base and ceiling of the layer
         :param h: Height at which to evaluate the temperature, pressure, density
-        :return:
+        :return: temperature, pressure, density, speed of sound
         """
+
         h_max = self.get_ceiling_height()
+        h0 = self.get_base_height()
 
         if h > h_max:
             raise ValueError
 
-        if h == self.get_base_height():
+        if h == h0:
             return self.get_base_values()
 
-        h0, T0, P0, D0 = self.get_base_values()
+        T0, P0, D0, a0, mu0 = self.get_base_values()
 
         L = (self.__T_top - T0) / (h_max - h0)
         C = -self.g0 / (L * self.R)  # To Simplify and shorten code we define the following expression for the exponent
@@ -182,5 +218,8 @@ def get_intermediate_values(self, h) -> list:
         P = P0 * (T / T0) ** C
         D = D0 * (T / T0) ** (C - 1)
 
-        return [h, T, P, D]
+        a = self.speed_of_sound(T)
+        mu = self.sutherland_viscosity(T)
+
+        return [T, P, D, a, mu]
 

setup.py

@@ -9,14 +9,14 @@ def readme():
 setup(
   name = 'isacalc',         # How you named your package folder (MyLib)
   packages = ['isacalc'],   # Chose the same as "name"
-  version = 'v1.1.0',      # Start with a small number and increase it with every change you make
+  version = 'v1.2.0',      # Start with a small number and increase it with every change you make
   license='GNU GPLv3',        # Chose a license from here: https://help.github.com/articles/licensing-a-repository
   description = 'Standard International Atmosphere Calculator',   # Give a short description about your library
   long_description=readme(),
   author = 'Luke de Waal',                   # Type in your name
   author_email = 'lr.de.waal.01@gmail.com',      # Type in your E-Mail
   url = 'https://github.com/LukeDeWaal/ISA_Calculator',   # Provide either the link to your github or to your website
-  download_url = 'https://github.com/LukeDeWaal/ISA_Calculator/archive/v1.1.0.tar.gz',    # I explain this later on
+  download_url = 'https://github.com/LukeDeWaal/ISA_Calculator/archive/v1.2.0.tar.gz',    # I explain this later on
   keywords = ['ISA','Aerospace','Aeronautical','Atmosphere'],   # Keywords that define your package best
   install_requires=[            # I get to this in a second
           'numpy',