Initial commit

.github/workflows/Flake8.yml

@@ -0,0 +1,30 @@
+name: Linting
+
+on:
+  push:
+    branches: 
+      - master
+
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v1
+    - name: Set up Python 3.6
+      uses: actions/setup-python@v1
+      with:
+        python-version: 3.6
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install -r requirements.txt
+        
+    - name: Lint with flake8
+      run: |
+        pip install flake8
+        # stop the build if there are Python syntax errors or undefined names
+        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
+        # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
+        flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics

.github/workflows/UnitTests.yml

@@ -0,0 +1,34 @@
+name: Tests
+
+on:
+  push:
+    branches: 
+      - master
+
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v1
+    - name: Set up Python 3.8
+      uses: actions/setup-python@v1
+      with:
+        python-version: 3.8
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install -r requirements.txt
+        pip install coverage
+        python setup.py build_ext --inplace
+                
+    - name: Test with unittest
+      run: |
+        export PYTHONPATH="$PYTHONPATH:$(pwd)/MAPLEAF"
+        python -m coverage run --source=./MAPLEAF -m unittest discover -v
+        
+    - name: Upload coverage results
+      run: |
+        export CODECOV_TOKEN="5b4f14f3-29a8-4f93-96e7-4723ecb7749a"
+        bash <(curl -s https://codecov.io/bash)

.gitignore

@@ -0,0 +1,79 @@
+# Ignore the testWrite file, it will have a new date everytime the unit tests are run
+testWrite.mapleaf
+testSimConfigWrite.mapleaf
+TempFigureDir
+
+#Ignore Pycache
+__pycache__
+
+#Remove Python intermediate files
+*.pyc
+/dist/
+/*.egg-info
+build/
+
+#Remove .vscode files
+*.vscode
+
+#Notcom output files
+for00[123456789].dat
+for01*
+
+#Simulation Logs
+*simulationLog*
+*forceEvaluationLog*
+*monteCarloLog*
+
+# Ignore pypy environment
+my-pypy-env
+
+# Ignore profiling outputs
+*.prof
+~lock*
+
+# Ignore code coverage output
+.coverage
+
+# Ignore convergence result files
+convergenceResult.csv
+adaptiveConvergenceResult.csv
+adaptiveConvergenceResult.xlsx
+convergenceResult.xlsx
+~*
+
+# Exclude coverage info
+htmlcov
+coverage.xml
+
+# Exclude class/package diagrams
+packages*.png
+packages*.dot
+packages*.pdf
+classes*.png
+classes*.dot
+classes*.pdf
+
+# Exclude regression testing plots
+test/regressionTesting/*.pdf
+test/regressionTesting/*.png
+test/regressionTesting/*.eps
+test/regressionTesting/*/*.pdf
+test/regressionTesting/*/*.png
+test/regressionTesting/*/*.eps
+test/regressionTesting/*/*/*.pdf
+test/regressionTesting/*/*/*.png
+test/regressionTesting/*/*/*.eps
+
+# Exlude temporary files generated with regression testing results in them
+*_newExpectedResultsRecorded*
+
+# Ignore Cython intermediate .c files & compiled files
+*Cython*.c
+*Cython*.cpp
+*Cython*.pyd
+*Cython*.so
+
+clean.sh
+
+# Documentation
+doc

CodingStyle.md

@@ -0,0 +1,37 @@
+## Document outlining coding practices and style guidelines
+
+## Naming
+### Variables:
+* mixedCase (from PEP8)
+* Use descriptive names x = bad, numEngines = good
+
+### Classes:
+* CapitalizedWords (from PEP8)
+
+### Modules:
+* Module names should match class name when a module contains a single class (most common case) (ex. Rocket.py, Vector.py)
+* Module name should be plural or a verb if it holds multiple classes / functions (ex. MeanWindModelling.py, RocketComponents.py)
+
+### Methods/Functions:
+* Use a leading underscore to indicate that a method (or variable) is only intended for internal use (private)
+
+### Underscores:
+* Use underscores when a clear separation is required in a name
+* Used most commonly in this repository to separate a common prefix from unique suffixes in names (ex. RigidBody_3DoF / RigidBody_6DoF)
+
+## Formatting
+* An auto-formatter is not being used yet. May switch to one in the future  
+* Try to wrap long lines (> 120 characters) . For long/verbose function calls/definitions, put each argument on its own line  
+* In math, whenever possible, use spacing to hint at the order of operations, Ex: `(x*x + y*y)` is preferred over `(x * x + y * y)`  
+* No blank lines immediately after control flow statements
+
+## Other
+### Unit Tests:
+* Modules should be in ./test/
+* Module names should begin with test_
+* Class names should begin with Test (and inherit from unittest.TestCase)
+* Test method names should begin with test_
+
+### Notes:
+* Try and keep code as modular as possible so that changes and improvements can be made easily. Create/specify concise interfaces. Use abstract base classes (abc.ABC) to specify interfaces
+* Use functions and classes to simplify code and make its intentions obvious

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2020 Henry Stoldt and other contributors
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

MAPLEAF/ENV/AtmosphereModelling.py

@@ -0,0 +1,177 @@
+''' Atmospheric models used by `MAPLEAF.ENV.Environment.Environment` defined and initialized here '''
+
+import abc
+from bisect import bisect
+from math import exp
+from typing import Sequence
+
+import numpy as np
+
+from MAPLEAF.Interpolation import linInterp
+from MAPLEAF.IO.SimDefinition import defaultConfigValues
+
+
+class AtmosphericModel(abc.ABC):
+    ''' Defines the interface of atmospheric models returned by atmosphericModelFactory '''
+
+    @abc.abstractmethod
+    def getAirProperties(self, ASLElevation: float, time: float) -> Sequence[float]:
+        ''' 
+            Should return an iterable containing: 
+                temp(K), 
+                static pressure (Pa), 
+                density (kg/m^3), 
+                dynamic viscosity (Pa*s),
+            in that order
+        '''
+        return
+
+def atmosphericModelFactory(atmosphericModel=None, envDictReader=None):
+    ''' 
+        Provide either an atmosphericModel name ('USStandardAtmosphere' is only option right now that doesn't require additional info,
+            or provide an envDictReader (`MAPLEAF.IO.SubDictReader.SubDictReader`)
+    '''
+    if atmosphericModel == None:
+        atmosphericModel = envDictReader.getString("AtmosphericPropertiesModel")
+
+    if atmosphericModel == "Constant":
+        if envDictReader == None:
+            raise ValueError("envDictReader required to initialize Constant atm properties model")
+
+        # Get values from ConstantAtmosphere subDictionary
+        envDictReader.simDefDictPathToReadFrom = "Environment.ConstantAtmosphere"
+
+        constTemp = envDictReader.getFloat("temp") + 273.15 # Convert to Kelvin (Expecting Celsius input)
+        constPressure = envDictReader.getFloat("pressure")
+        constDensity = envDictReader.getFloat("density")
+        constViscosity = envDictReader.getFloat("viscosity")
+
+        # Return to reading from Environment for any subsequent parsing
+        envDictReader.simDefDictPathToReadFrom = "Environment"
+
+        return ConstantAtmosphericProperties(constTemp, constPressure, constDensity, constViscosity)
+
+    elif atmosphericModel == "TabulatedAtmosphere":
+        try:
+            tableFilePath = envDictReader.getString("TabulatedAtmosphere.filePath")
+        except AttributeError:
+            tableFilePath = defaultConfigValues["Environment.TabulatedAtmosphere.filePath"]
+        return TabulatedAtmosphere(tableFilePath)
+
+    elif atmosphericModel == "USStandardAtmosphere":
+        return USStandardAtmosphere()
+
+    else:
+        raise ValueError("Atmospheric model: {} not implemented, try using 'USStandardAtmosphere'".format(atmosphericModel))
+
+class ConstantAtmosphericProperties(AtmosphericModel):
+    def __init__(self, temp, pressure, density, viscosity):
+        self.airProperties = [ temp, pressure, density, viscosity ]
+
+    def getAirProperties(self, _, _2=None):
+        return self.airProperties
+
+class TabulatedAtmosphere(AtmosphericModel):
+    '''
+        Provides linearly-interpolated atmospheric properties from a table in a file
+        Table columns are expected are:
+            h(m ASL)	T(K)	P(Pa)	rho(kg/m^3)	    mu(10^-5 Pa*s)
+    '''
+    def __init__(self, filePath):
+        inputTable = np.loadtxt(filePath, skiprows=1)
+
+        # Convert viscosity to Pa-s
+        inputTable[:, 4] *= 1E-5 
+
+        # Set up keys, values for interpolation of all properties at once
+        self.keys = inputTable[:, 0] # Altitude (m ASL) is interpolation key
+        self.values = inputTable[:, 1:] # Contains: T(K), P(Pa), rho(kg/m^3), dynamicViscosity(Pa*s), all to be interpolated over
+
+    def getAirProperties(self, ASLElevation, _=None):
+        return linInterp(self.keys, self.values, ASLElevation)
+
+class USStandardAtmosphere(AtmosphericModel):
+    '''
+        Provides atmospheric properties calculated according to the model here: 
+        https://nebula.wsimg.com/ab321c1edd4fa69eaa94b5e8e769b113?AccessKeyId=AF1D67CEBF3A194F66A3&disposition=0&alloworigin=1
+
+    '''
+    T0 = 288.15 # K (15 Celsius) @ ASL == 0
+    p0 = 101325 # Pa @ ASL == 0
+    M = 28.9644 # Molecular weight of air
+    R = 8.31432 # J/molK
+    earthRadius = 6356766 # m - used to convert to geopotential altitude
+    G = 9.80665 # m/s^2 (sea level, 45 degree latitude)
+
+    baseHeights = [ -2000, 11000, 20000, 32000, 47000, 51000, 71000, 84852 ] # m (Geopotential)
+    dt_dh = [ -6.5e-3, 0, 1.0e-3, 2.8e-3, 0, -2.8e-3, -2.0e-3, 0 ] # K/m
+    baseTemps = [] # Gets filled out in self.__init__
+    basePressures = [] # Gets filled out in self.__init__
+
+    def __init__(self):
+        # Pre-Calculate base temperatures for each interval
+        dt_dh = self.dt_dh[0]
+        baseTemp1 = self.T0 + dt_dh*self.baseHeights[0]
+        self.baseTemps.append(baseTemp1)
+
+        Pb_over_P0 = ( (self.T0 + dt_dh*self.baseHeights[0]) / \
+            self.T0) ** (-self.G*self.M/(self.R * dt_dh * 1000))
+        self.basePressures.append(Pb_over_P0 * self.p0)
+
+        # Loop over each layer, calculate temp & pressure at base
+        for i in range(1, len(self.baseHeights)):
+            # Calculate temperatures at base of each layer
+            Tb, Pb, dt_dh = self.baseTemps[-1], self.basePressures[i-1], self.dt_dh[i-1]
+            dh = float(self.baseHeights[i] - self.baseHeights[i-1])
+
+            nextTemp = Tb + dt_dh*dh
+
+            # Calculate pressures at the base of each layer
+            if dt_dh == 0.0:
+                P_over_Pb = exp(-self.G*self.M*dh / (self.R * Tb * 1000))
+                nextPressure = P_over_Pb * Pb
+
+            else:
+                body = (Tb + dt_dh*dh) / Tb
+                exponent = -self.G*self.M / (self.R*dt_dh*1000)
+                nextPressure = body**exponent * Pb
+
+            self.baseTemps.append(nextTemp)
+            self.basePressures.append(nextPressure)
+
+    def getAirProperties(self, ASLElevation, time):
+
+        # Calculate geopotential altitude (H)
+        H = (self.earthRadius * ASLElevation) / (self.earthRadius + ASLElevation)
+
+        # Set density to zero above 86 km
+        if H >= 86000:
+            H = 86000
+
+        # Figure out which interval we're in 
+        baseIndex = bisect(self.baseHeights, H) - 1
+        altitudeAboveBase = H - self.baseHeights[baseIndex] # Geopotential altitude above base
+        dt_dh = self.dt_dh[baseIndex]
+
+        # Calc temp
+        Tb = self.baseTemps[baseIndex]
+        temp = Tb + dt_dh * altitudeAboveBase
+
+        # Calc pressure
+        Pb = self.basePressures[baseIndex]
+        if dt_dh == 0:
+            P_over_Pb = exp(-self.G*self.M * altitudeAboveBase / (self.R * Tb * 1000))
+            pressure = P_over_Pb * Pb
+        else:
+            body = (Tb + dt_dh*altitudeAboveBase) / Tb
+            exponent = -self.G*self.M / (self.R*dt_dh*1000)
+            pressure = body**exponent * Pb
+
+        # Calculate density and viscosity from temperature and pressure   
+        rho = self.M * pressure / (self.R * temp * 1000) # Ideal gas law
+        viscosity = 1.457e-6 * temp**(1.5) / (temp + 110.4) # Sutherland's law
+
+        if H >= 86000:
+            rho = 0.0
+
+        return [ temp, pressure, rho, viscosity ]