Module `MAPLEAF.ENV`

Environmental modelling: main class is Environment. Environment wraps the atmospheric models, gravity models, mean wind models, and turbulence models.

Expand source code

'''
Environmental modelling: main class is `MAPLEAF.ENV.Environment`.
`MAPLEAF.ENV.Environment` wraps the atmospheric models, gravity models, mean wind models, and turbulence models.

.. image:: https://www.goodfreephotos.com/albums/astrophotography/sunrise-over-the-earth.jpg
'''
# Make the classes in all submodules importable directly from MAPLEAF.Rocket
from .launchRail import *
from .MeanWindModelling import *
from .EarthModelling import *
from .AtmosphereModelling import *
from .TurbulenceModelling import *
from .environment import *

subModules = [ environment, launchRail, MeanWindModelling, EarthModelling, AtmosphereModelling, TurbulenceModelling ]

__all__ = []

for subModule in subModules:
    __all__ += subModule.__all__

Sub-modules

MAPLEAF.ENV.AtmosphereModelling: These classes model the change of air properties (Pressure, Density, etc… with altitude)
MAPLEAF.ENV.EarthModelling: These classes model the earth's gravity and perform some coordinate transformations between the global inertial frame and the ENU/launch tower frame
MAPLEAF.ENV.MeanWindModelling: Modeling of the mean / average component of wind velocity
MAPLEAF.ENV.TurbulenceModelling: Modeling of turbulent, fluctuating component of wind velocity
MAPLEAF.ENV.environment: Main wrapper and data classes that tie together all of the environmental models and are queried by instances of Rocket
MAPLEAF.ENV.launchRail: Modeling of the effects of a Launch Rail on rocket launches

Functions

def atmosphericModelFactory(envDictReader=None) ‑> MAPLEAF.ENV.AtmosphereModelling.AtmosphericModel

Provide either an atmosphericModel name ('USStandardAtmosphere' is only option right now that doesn't require additional info, or provide an envDictReader (SubDictReader)

Expand source code

def atmosphericModelFactory(envDictReader=None) -> AtmosphericModel:
    ''' 
        Provide either an atmosphericModel name ('USStandardAtmosphere' is only option right now that doesn't require additional info,
            or provide an envDictReader (`MAPLEAF.IO.SubDictReader`)
    '''
    if envDictReader == None:
        modelType = defaultConfigValues["Environment.AtmosphericPropertiesModel"]
    else:
        modelType = envDictReader.getString("AtmosphericPropertiesModel")

    if modelType == "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 ConstantAtmosphere(constTemp, constPressure, constDensity, constViscosity)
    
    elif modelType == "TabulatedAtmosphere":
        try:
            tableFilePath = envDictReader.getString("TabulatedAtmosphere.filePath")
        except AttributeError:
            tableFilePath = defaultConfigValues["Environment.TabulatedAtmosphere.filePath"]

        tableFilePath = getAbsoluteFilePath(tableFilePath)
        return TabulatedAtmosphere(tableFilePath)

    elif modelType == "USStandardAtmosphere":
        return USStandardAtmosphere()

    else:
        raise ValueError("Atmospheric model: {} not implemented, try using 'USStandardAtmosphere'".format(modelType))

def earthModelFactory(envDictReader=None) ‑> MAPLEAF.ENV.EarthModelling.EarthModel

Provide an envDictReader (SubDictReader). If none is provided, returns a FlatEarth model

Expand source code

def earthModelFactory(envDictReader=None) -> EarthModel:
    ''' Provide an envDictReader (`MAPLEAF.IO.SubDictReader`). If none is provided, returns a FlatEarth model '''
    if envDictReader == None:
        return FlatEarth()

    earthModel = envDictReader.getString("EarthModel")

    if earthModel == "Flat":
        return FlatEarth()
    elif earthModel == "Round":
        return SphericalEarth()
    elif earthModel == "WGS84":
        return WGS84()
    elif earthModel == "None":
        return NoEarth()
    else:
        raise NotImplementedError("Earth model: {} not found. Try 'Flat', 'Round' or 'WGS84'".format(earthModel))

def meanWindModelFactory(envReader=None, silent=False) ‑> MAPLEAF.ENV.MeanWindModelling.MeanWindModel

Instantiates a mean wind model

Expand source code

def meanWindModelFactory(envReader=None, silent=False) -> MeanWindModel:
    ''' Instantiates a mean wind model '''

    if envReader == None:
        # Return default zero-wind model
        constWind = Vector(defaultConfigValues["Environment.ConstantMeanWind.velocity"])
        return ConstantWind(constWind)

    meanWindModel = None
    meanWindModelType = envReader.getString("MeanWindModel")

    if meanWindModelType == "Constant":
        meanGroundWind = envReader.getVector("ConstantMeanWind.velocity")
        if not silent:
            print("Constant ground wind: {:1.2f} m/s".format(meanGroundWind))

        meanWindModel = ConstantWind(meanGroundWind)

    elif meanWindModelType in [ "SampledGroundWindData", "Hellman" ]:
        def getLocationSampledGroundWindVel():
            #### Get list of, and weights of, locations to sample wind rose from
            GroundWindLocationsSampled = envReader.getString("SampledGroundWindData.locationsToSample").split()
            locationsSampled = []
            locationWeights = []
            # Parse into locations (1st, 3rd, 5th... values) and weights (2nd, 4th, 6th... values)
            for i in range(0, len(GroundWindLocationsSampled), 2):
                locationsSampled.append(GroundWindLocationsSampled[i])
                locationWeights.append(float(GroundWindLocationsSampled[i+1]))
            launchMonth = envReader.getString("SampledGroundWindData.launchMonth")

            # Sample wind rose(s)
            sampler = WindRoseDataSampler(silent=silent)
            meanGroundWind = sampler.sampleWindRoses(locationsSampled, locationWeights, launchMonth)

            # Output choices parsed from input file and resulting wind
            if not silent:
                if envReader.simDefinition.monteCarloLogger != None:
                    envReader.simDefinition.monteCarloLogger.log("Sampling ground winds from: {}".format(locationsSampled))
                    envReader.simDefinition.monteCarloLogger.log("Sampling weights for each location: {}".format(locationWeights))
                    envReader.simDefinition.monteCarloLogger.log("Sampling wind distribution from month of: {}".format(launchMonth))
                    envReader.simDefinition.monteCarloLogger.log("Sampled mean ground wind: {:1.2f} m/s".format(meanGroundWind))
                else:
                    print("Sampling ground winds from: {}".format(locationsSampled))
                    print("Sampling weights for each location: {}".format(locationWeights))
                    print("Sampling wind distribution from month of: {}".format(launchMonth))
                    print("Sampled mean ground wind: {:1.2f} m/s".format(meanGroundWind))

            return meanGroundWind

        ### Create and return appropriate wind vs altitude model ###
        if meanWindModelType == "SampledGroundWindData":
            meanGroundWind = getLocationSampledGroundWindVel()
            if not silent:
                print("Wind is not a function of altitude")
            meanWindModel = ConstantWind(meanGroundWind)

        elif meanWindModelType == "Hellman":
            groundWindModel = envReader.getString("Hellman.groundWindModel")

            # Get ground wind
            if groundWindModel == "Constant":
                meanGroundWind = envReader.getVector("ConstantMeanWind.velocity")
            elif groundWindModel == "SampledGroundWindData":
                meanGroundWind = getLocationSampledGroundWindVel()

            HellmanAlphaCoeff = envReader.getFloat("Hellman.alphaCoeff")
            HellmanAltitudeLimit = envReader.getFloat("Hellman.altitudeLimit")

            if not silent:
                print("Constant ground wind: {:1.2f} m/s".format(meanGroundWind))
                print("Wind vs. altitude governed by Hellman law: v2 = v1*(z2/10)^a where a = {}".format(HellmanAlphaCoeff))
                print("Hellman law is taken to apply up to an altitude of {} m AGL".format(HellmanAltitudeLimit))

            meanWindModel = Hellman(meanGroundWind, HellmanAltitudeLimit, HellmanAlphaCoeff)
    
    elif meanWindModelType == "CustomWindProfile":
        meanWindProfileFilePath = envReader.getString("CustomWindProfile.filePath")
        meanWindModel = InterpolatedWind(windFilePath=meanWindProfileFilePath)

    elif meanWindModelType == "SampledRadioSondeData":
        # Get locations and location weights
        locationsAndWeights = envReader.getString("SampledRadioSondeData.locationsToSample").split()
        locations = []
        weights = []
        for i in range(0, len(locationsAndWeights), 2):
            locations.append(locationsAndWeights[i])
            weights.append(float(locationsAndWeights[i+1]))
            
        locationASLAltitudes = [ float(x) for x in envReader.getString("SampledRadioSondeData.locationASLAltitudes").split() ]

        # Get launch month
        launchMonth = envReader.getString("SampledRadioSondeData.launchMonth")

        # Get random seed (if provided)
        radioSondeRandomSeed = envReader.tryGetString("SampledRadioSondeData.randomSeed")

        # Select and parse radio sonde profile
        sampler = RadioSondeDataSampler(silent=silent)
        altitudes, windVectors = sampler.getRadioSondeWindProfile(locations, weights, locationASLAltitudes, launchMonth, radioSondeRandomSeed)
        
        meanWindModel = InterpolatedWind(windAltitudes=altitudes, winds=windVectors)

    else:
        raise ValueError("Unknown MeanWindModel: {}. Please see SimDefinitionTemplate.txt for available options.".format(meanWindModelType))

    return meanWindModel

def turbulenceModelFactory(envReader, silent=False)

Reads data from simDefinition, initializes and returns the appropriate TurbulenceModel

Expand source code

def turbulenceModelFactory(envReader, silent=False):
    ''' Reads data from simDefinition, initializes and returns the appropriate TurbulenceModel '''
    turbModelType = envReader.getString("Environment.TurbulenceModel")
    
    if not silent:
        print("Turbulence/Gust Model: {}".format(turbModelType))

    if turbModelType== "None":
        turbulenceModel = NoTurb()
        
    elif "PinkNoise" in turbModelType:
        measuredPinkNoiseStdDev = 2.26 # For a 2-pole PinkNoiseGenerator - re-measure if the number of poles is changed

        turbulenceIntensity = envReader.tryGetInt("PinkNoiseModel.turbulenceIntensity")
        if turbulenceIntensity != None:
            turbulenceIntensity /= (measuredPinkNoiseStdDev * 100)

        velocityStDev = envReader.tryGetInt("PinkNoiseModel.velocityStdDeviation")
        if velocityStDev != None:
            velocityStDev /= measuredPinkNoiseStdDev

        pinkNoiseRandomSeed1 = envReader.tryGetInt("PinkNoiseModel.randomSeed1")
        pinkNoiseRandomSeed2 = envReader.tryGetInt("PinkNoiseModel.randomSeed2")
        pinkNoiseRandomSeed3 = envReader.tryGetInt("PinkNoiseModel.randomSeed3")
            
        if turbModelType == "PinkNoise1D":
            turbulenceModel = PinkNoise1D(turbulenceIntensity, velocityStDev, pinkNoiseRandomSeed1)
            if not silent:
                print("Random seed 1: {}".format(turbulenceModel.png1.seed))

        elif turbModelType == "PinkNoise2D":
            turbulenceModel = PinkNoise2D(turbulenceIntensity, velocityStDev, pinkNoiseRandomSeed1, pinkNoiseRandomSeed2)
            if not silent:
                print("Random seed 1: {}".format(turbulenceModel.png1.seed))
                print("Random seed 2: {}".format(turbulenceModel.png2.seed))

        elif turbModelType == "PinkNoise3D":
            turbulenceModel = PinkNoise3D(turbulenceIntensity, velocityStDev, pinkNoiseRandomSeed1, pinkNoiseRandomSeed2, pinkNoiseRandomSeed3)
            if not silent:
                print("Random seed 1: {}".format(turbulenceModel.png1.seed))
                print("Random seed 2: {}".format(turbulenceModel.png2.seed))
                print("Random seed 2: {}".format(turbulenceModel.png3.seed))
            
    elif turbModelType == "customSineGust":
        GustStartAltitude = envReader.getFloat("CustomSineGust.startAltitude")
        GustMagnitude = envReader.getFloat("CustomSineGust.magnitude")
        GustSineBlendDistance = envReader.getFloat("CustomSineGust.sineBlendDistance")
        GustLayerThickness = envReader.getFloat("CustomSineGust.thickness")
        GustDirection = envReader.getVector("CustomSineGust.direction").normalize()

        turbulenceModel = CustomSineGust(GustStartAltitude, GustSineBlendDistance, GustLayerThickness, GustMagnitude, GustDirection)

        if not silent:
            print("Gust Altitude: {}".format(GustStartAltitude))
            print("Gust Vector: {}".format(GustMagnitude * GustDirection))

    else:
        raise ValueError("Turbulence Model {} not found. Please choose from: None, PinkNoise1/2/3D or customSineGust")

    return turbulenceModel

Classes

class ConstantAtmosphere (temp, pressure, density, viscosity)

Interface for all atmosphere models

Expand source code

class ConstantAtmosphere(AtmosphericModel):
    def __init__(self, temp, pressure, density, viscosity):
        self.airProperties = [ temp, pressure, density, viscosity ]

    def getAirProperties(self, _, _2=None):
        return self.airProperties

Ancestors

MAPLEAF.ENV.AtmosphereModelling.AtmosphericModel
abc.ABC

Methods

def getAirProperties(self, _)

Return an iterable containing: temp(K), static pressure (Pa), density (kg/m^3), dynamic viscosity (Pa*s), in that order

Expand source code

def getAirProperties(self, _, _2=None):
    return self.airProperties

class ConstantWind (wind)

Defines a constant wind speed at all altitudes

Expand source code

class ConstantWind(MeanWindModel):
    ''' Defines a constant wind speed at all altitudes  '''
    
    def __init__(self, wind):
        self.wind = wind

    def getMeanWind(self, AGLAltitude):
        return self.wind

Ancestors

MAPLEAF.ENV.MeanWindModelling.MeanWindModel
abc.ABC

Methods

def getMeanWind(self, AGLAltitude)

Expand source code

def getMeanWind(self, AGLAltitude):
    return self.wind

class Environment (simDefinition=None, silent=False)

Class wraps Wind Models, Atmospheric properties models, and earth/gravity models, presenting a single interface for communication with flight vehicles

Sets up the Wind, Atmospheric, and Earth models requested in the Sim Definition file

Expand source code

class Environment():
    '''
        Class wraps Wind Models, Atmospheric properties models, and earth/gravity models, presenting a 
        single interface for communication with flight vehicles
    '''
    
    def __init__(self, simDefinition=None, silent=False):
        ''' Sets up the Wind, Atmospheric, and Earth models requested in the Sim Definition file '''
        self.launchRail = None

        if simDefinition != None:
            # Whenever we're running a real simulation, should always end up here
            envDictReader = SubDictReader("Environment", simDefinition)

            self.meanWindModel = meanWindModelFactory(envDictReader, silent=silent)
            self.turbulenceModel = turbulenceModelFactory(envDictReader, silent=silent)
            self.atmosphericModel = atmosphericModelFactory(envDictReader=envDictReader)
            self.earthModel = earthModelFactory(envDictReader)

            self.launchSiteElevation = envDictReader.tryGetFloat("LaunchSite.elevation")
            self.launchSiteLatitude = envDictReader.tryGetFloat("LaunchSite.latitude")
            self.launchSiteLongitude = envDictReader.tryGetFloat("LaunchSite.longitude")

            # Check if being launched from a launch rail
            launchRailLength = envDictReader.getFloat("LaunchSite.railLength")
            
            if launchRailLength > 0:
                # Initialize a launch rail, aligned with the rocket's initial direction
                initialRocketPosition_towerFrame = envDictReader.getVector("Rocket.position")

                # Check whether precise initial orientation has been specified
                rotationAxis = envDictReader.tryGetVector("Rocket.rotationAxis", defaultValue=None)
                if rotationAxis != None:
                    rotationAngle = math.radians(envDictReader.getFloat("Rocket.rotationAngle"))
                    initOrientation = Quaternion(rotationAxis, rotationAngle)
                else:
                    # Calculate initial orientation quaternion in launch tower frame
                    rotationAxis = envDictReader.tryGetVector("Rocket.rotationAxis", defaultValue=None)
                    if rotationAxis != None:
                        rotationAngle = math.radians(self.rocketDictReader.getFloat("Rocket.rotationAngle"))
                        initOrientation = Quaternion(rotationAxis, rotationAngle)
                    else:
                        # Calculate initial orientation quaternion in launch tower frame
                        initialDirection = envDictReader.getVector("Rocket.initialDirection").normalize()
                        angleFromVertical = Vector(0,0,1).angle(initialDirection)
                        rotationAxis = Vector(0,0,1).crossProduct(initialDirection)
                        initOrientation = Quaternion(rotationAxis, angleFromVertical)

                    # TODO: Get from rocket, or calculate the same way - so that it works with rotationAxis + Angle
                    initialDirection = envDictReader.getVector("Rocket.initialDirection").normalize()
                    angleFromVertical = Vector(0,0,1).angle(initialDirection)
                    rotationAxis = Vector(0,0,1).crossProduct(initialDirection)
                    initOrientation = Quaternion(rotationAxis, angleFromVertical)

                initPosition_seaLevelENUFrame = initialRocketPosition_towerFrame + Vector(0, 0, self.launchSiteElevation)
                launchTowerState_local = RigidBodyState(position=initPosition_seaLevelENUFrame, orientation=initOrientation)
                launchTowerState_global = self.earthModel.convertIntoGlobalFrame(launchTowerState_local, self.launchSiteLatitude, self.launchSiteLongitude)
                towerDirection_global = launchTowerState_global.orientation.rotate(Vector(0, 0, 1))                
                self.launchRail = LaunchRail(launchTowerState_global.position, towerDirection_global, launchRailLength, earthRotationRate=self.earthModel.rotationRate)
                
        else:
            # Construct default environment from default parameters when no sim definition is passed in
            # Need additional default value here in case a SimDefinition is not passed in (otherwise SimDefinition takes care of default values)
            self.meanWindModel = meanWindModelFactory()
            self.turbulenceModel = None
            self.atmosphericModel = atmosphericModelFactory()
            self.earthModel = earthModelFactory()

            self.launchSiteElevation = float(defaultConfigValues["Environment.LaunchSite.elevation"])
            self.launchSiteLatitude = float(defaultConfigValues["Environment.LaunchSite.latitude"])
            self.launchSiteLongitude = float(defaultConfigValues["Environment.LaunchSite.longitude"])

    def convertInitialStateToGlobalFrame(self, initialStateInLaunchTowerFrame):
        ''' 
            Used to convert the rocket's initial kinematic state from the launch tower frame (in which it is defined) and 
                into the global inertia frame

            Takes a rigid body state defined in the launch tower frame and converts it to the global frame 
            For a flat earth model, this just adjusts for the altitude of the launch site,
            but for a rotating earth model, this modifies every part of the rigid body state:
                1. Position is redefined relative to the center of the earth (Acoording to lat/lon)  
                2. The velocity of earth's rotation is added to the velocity  
                3. The orientation is redefined relative to the global frame  
                4. The angular velocity of the earth is added to the rocket's initial state  

            .. note:: If a launch rail is being used, set the launch-tower-frame angular velocity to zero before doing the conversion
        '''
        # In all cases, first redefine the present state relative to sea level
        initialStateInLaunchTowerFrame.position += Vector(0,0,self.launchSiteElevation)

        # Set launch-tower-frame angular velocity to zero if using a launch rail
        if self.launchRail != None:
            initialStateInLaunchTowerFrame.angularVelocity = AngularVelocity(0,0,0)

        # Call the current earth model's conversion function
        return self.earthModel.convertIntoGlobalFrame(initialStateInLaunchTowerFrame, self.launchSiteLatitude, self.launchSiteLongitude)

    def convertStateToENUFrame(self, globalFrameState: Union[RigidBodyState, RigidBodyState_3DoF]) -> Union[RigidBodyState, RigidBodyState_3DoF]:
        altitude = self.earthModel.getAltitude(*globalFrameState.position)
        position = Vector(0, 0, altitude) # Frame moves with the aircraft so x/y are always zero
        
        inertialToENURotation = self.earthModel.getInertialToENUFrameRotation(*globalFrameState.position)
        ENUToGlobalRotation = inertialToENURotation.conjugate()
        velocity = ENUToGlobalRotation.rotate(globalFrameState.velocity)

        try:
            orientation = ENUToGlobalRotation * globalFrameState.orientation
            angVel = ENUToGlobalRotation.rotate(globalFrameState.angularVelocity)
            return RigidBodyState(position, velocity, orientation, angVel)
        except:
            return RigidBodyState_3DoF(position, velocity)

    #### Get all air/atmospheric properties ####
    def getAirProperties(self, position: Vector, time=None) -> EnvironmentalConditions:
        ''' Pass in a vector representing the aircraft's position in the global inertial frame of reference '''
        # ASL Altitude is calculated differently depending on whether the earth is modelled as flat, round, or ellipsoidal
        ASLAltitude = self.earthModel.getAltitude(*position) 
        # Ground level is remain constant as that of the launch site
        AGLAltitude = ASLAltitude - self.launchSiteElevation

        # Get wind based on AGL, in North-East-Up frame
            # TODO: This is only desirable for low-altitude flights, should switch to ASL altitude once out of the atmospheric boundary layer
        meanWind, turbWind = self.getWind(AGLAltitude, time)
        
        # Add earth's rotation speed to the mean wind
        # For No/Flat Earth models, the rotationRate will be zero and this will have no effect
        distanceFromRotationAxis = math.sqrt(position.X**2 + position.Y**2)
        earthRotationSpeed = self.earthModel.rotationRate * distanceFromRotationAxis
        # Earth's surface rotates towards the east, so add the rotation velocity to the East component (X) of the local ENU frame
        meanWind = meanWind + Vector(earthRotationSpeed,0,0)

        # Wind is calculated in the North-East-Up frame, and needs to be rotated into the global inertial frame
        # For No/Flat Earth models, the rotation will return a zero-rotation Quaternion, so this has no effect
        orientationOfENUFrameInGlobalFrame = self.earthModel.getInertialToENUFrameRotation(*position)
        meanWind = orientationOfENUFrameInGlobalFrame.rotate(meanWind)
        turbWind = orientationOfENUFrameInGlobalFrame.rotate(turbWind)

        # Now calculate total wind in global frame
        totalWind = meanWind + turbWind

        # Add air properties based on ASL, return the result
        return EnvironmentalConditions(
            ASLAltitude,
            *self.atmosphericModel.getAirProperties(ASLAltitude, time), # Asterisk unpacks the returned values
            totalWind, 
            meanWind, 
            turbWind
        )

    #### Wind properties ####
    def getWind(self, AGLElevation, time=None):
        ''' Returns wind in the North-East-Up (Y-X-Z) frame '''
        meanWindVel = self.meanWindModel.getMeanWind(AGLElevation)
        
        if self.turbulenceModel != None:
            turbWindVel = self.turbulenceModel.getTurbVelocity(AGLElevation, meanWindVel, time)
        else:
            turbWindVel = Vector(0,0,0)

        return meanWindVel, turbWindVel

    #### Wrappers for launch rail functions ####    
    def applyLaunchRailMotionConstraints(self, state, time):
        if self.launchRail == None:
            # Majority of the time we'll be off the rail and this will run
            return state

        onLaunchRail, adjustedState = self.launchRail.applyLaunchTowerMotionConstraint(state, time)
        
        # Delete rail if we've left it
        if not onLaunchRail:
            print("Launch rail cleared")
            self.launchRail = None

        return adjustedState

    def applyLaunchTowerForce(self, state, time, unadjustedForce):
        if self.launchRail == None:
            # Majority of the time we'll be off the rail and this will run
            return unadjustedForce
        
        return self.launchRail.applyLaunchTowerForce(state, time, unadjustedForce)

    #### Gravity ####
    def getGravityForce(self, inertia, state) -> ForceMomentSystem:
        '''
            Inputs:
                inertia: (`MAPLEAF.Motion.Inertia`)
                state:   (`MAPLEAF.Motion.RigidBodyState`/`MAPLEAF.Motion.RigidBodyState_3DoF`)

            Returns:
                gravityForce: (ForceMomentSystem) defined in the rocket's local frame
        '''
        # Get gravity force in the global frame
        gravityForce = self.earthModel.getGravityForce(inertia, state)
        
        try:
            # Convert to local frame if in 6DoF simulation
            gravityForce = state.orientation.conjugate().rotate(gravityForce) # rotate into local frame when in 6DoF mode
        except AttributeError:
            pass # Don't do anything in 3DoF mode (No local frame exists)

        return ForceMomentSystem(gravityForce, inertia.CG)

Methods

def applyLaunchRailMotionConstraints(self, state, time)

Expand source code

def applyLaunchRailMotionConstraints(self, state, time):
    if self.launchRail == None:
        # Majority of the time we'll be off the rail and this will run
        return state

    onLaunchRail, adjustedState = self.launchRail.applyLaunchTowerMotionConstraint(state, time)
    
    # Delete rail if we've left it
    if not onLaunchRail:
        print("Launch rail cleared")
        self.launchRail = None

    return adjustedState

def applyLaunchTowerForce(self, state, time, unadjustedForce)

Expand source code

def applyLaunchTowerForce(self, state, time, unadjustedForce):
    if self.launchRail == None:
        # Majority of the time we'll be off the rail and this will run
        return unadjustedForce
    
    return self.launchRail.applyLaunchTowerForce(state, time, unadjustedForce)

def convertInitialStateToGlobalFrame(self, initialStateInLaunchTowerFrame)

Used to convert the rocket's initial kinematic state from the launch tower frame (in which it is defined) and into the global inertia frame

Takes a rigid body state defined in the launch tower frame and converts it to the global frame For a flat earth model, this just adjusts for the altitude of the launch site, but for a rotating earth model, this modifies every part of the rigid body state: 1. Position is redefined relative to the center of the earth (Acoording to lat/lon)
2. The velocity of earth's rotation is added to the velocity
3. The orientation is redefined relative to the global frame
4. The angular velocity of the earth is added to the rocket's initial state

Note: If a launch rail is being used, set the launch-tower-frame angular velocity to zero before doing the conversion

Expand source code

def convertInitialStateToGlobalFrame(self, initialStateInLaunchTowerFrame):
    ''' 
        Used to convert the rocket's initial kinematic state from the launch tower frame (in which it is defined) and 
            into the global inertia frame

        Takes a rigid body state defined in the launch tower frame and converts it to the global frame 
        For a flat earth model, this just adjusts for the altitude of the launch site,
        but for a rotating earth model, this modifies every part of the rigid body state:
            1. Position is redefined relative to the center of the earth (Acoording to lat/lon)  
            2. The velocity of earth's rotation is added to the velocity  
            3. The orientation is redefined relative to the global frame  
            4. The angular velocity of the earth is added to the rocket's initial state  

        .. note:: If a launch rail is being used, set the launch-tower-frame angular velocity to zero before doing the conversion
    '''
    # In all cases, first redefine the present state relative to sea level
    initialStateInLaunchTowerFrame.position += Vector(0,0,self.launchSiteElevation)

    # Set launch-tower-frame angular velocity to zero if using a launch rail
    if self.launchRail != None:
        initialStateInLaunchTowerFrame.angularVelocity = AngularVelocity(0,0,0)

    # Call the current earth model's conversion function
    return self.earthModel.convertIntoGlobalFrame(initialStateInLaunchTowerFrame, self.launchSiteLatitude, self.launchSiteLongitude)

def convertStateToENUFrame(self, globalFrameState: Union[RigidBodyState, RigidBodyState_3DoF]) ‑> Union[RigidBodyState, RigidBodyState_3DoF]

Expand source code

def convertStateToENUFrame(self, globalFrameState: Union[RigidBodyState, RigidBodyState_3DoF]) -> Union[RigidBodyState, RigidBodyState_3DoF]:
    altitude = self.earthModel.getAltitude(*globalFrameState.position)
    position = Vector(0, 0, altitude) # Frame moves with the aircraft so x/y are always zero
    
    inertialToENURotation = self.earthModel.getInertialToENUFrameRotation(*globalFrameState.position)
    ENUToGlobalRotation = inertialToENURotation.conjugate()
    velocity = ENUToGlobalRotation.rotate(globalFrameState.velocity)

    try:
        orientation = ENUToGlobalRotation * globalFrameState.orientation
        angVel = ENUToGlobalRotation.rotate(globalFrameState.angularVelocity)
        return RigidBodyState(position, velocity, orientation, angVel)
    except:
        return RigidBodyState_3DoF(position, velocity)

def getAirProperties(self, position: Vector, time=None) ‑> EnvironmentalConditions

Pass in a vector representing the aircraft's position in the global inertial frame of reference

Expand source code

def getAirProperties(self, position: Vector, time=None) -> EnvironmentalConditions:
    ''' Pass in a vector representing the aircraft's position in the global inertial frame of reference '''
    # ASL Altitude is calculated differently depending on whether the earth is modelled as flat, round, or ellipsoidal
    ASLAltitude = self.earthModel.getAltitude(*position) 
    # Ground level is remain constant as that of the launch site
    AGLAltitude = ASLAltitude - self.launchSiteElevation

    # Get wind based on AGL, in North-East-Up frame
        # TODO: This is only desirable for low-altitude flights, should switch to ASL altitude once out of the atmospheric boundary layer
    meanWind, turbWind = self.getWind(AGLAltitude, time)
    
    # Add earth's rotation speed to the mean wind
    # For No/Flat Earth models, the rotationRate will be zero and this will have no effect
    distanceFromRotationAxis = math.sqrt(position.X**2 + position.Y**2)
    earthRotationSpeed = self.earthModel.rotationRate * distanceFromRotationAxis
    # Earth's surface rotates towards the east, so add the rotation velocity to the East component (X) of the local ENU frame
    meanWind = meanWind + Vector(earthRotationSpeed,0,0)

    # Wind is calculated in the North-East-Up frame, and needs to be rotated into the global inertial frame
    # For No/Flat Earth models, the rotation will return a zero-rotation Quaternion, so this has no effect
    orientationOfENUFrameInGlobalFrame = self.earthModel.getInertialToENUFrameRotation(*position)
    meanWind = orientationOfENUFrameInGlobalFrame.rotate(meanWind)
    turbWind = orientationOfENUFrameInGlobalFrame.rotate(turbWind)

    # Now calculate total wind in global frame
    totalWind = meanWind + turbWind

    # Add air properties based on ASL, return the result
    return EnvironmentalConditions(
        ASLAltitude,
        *self.atmosphericModel.getAirProperties(ASLAltitude, time), # Asterisk unpacks the returned values
        totalWind, 
        meanWind, 
        turbWind
    )

def getGravityForce(self, inertia, state) ‑> ForceMomentSystem

Inputs

inertia: (Inertia) state: (RigidBodyState/RigidBodyState_3DoF)

Returns

gravityForce: (ForceMomentSystem) defined in the rocket's local frame

Expand source code

def getGravityForce(self, inertia, state) -> ForceMomentSystem:
    '''
        Inputs:
            inertia: (`MAPLEAF.Motion.Inertia`)
            state:   (`MAPLEAF.Motion.RigidBodyState`/`MAPLEAF.Motion.RigidBodyState_3DoF`)

        Returns:
            gravityForce: (ForceMomentSystem) defined in the rocket's local frame
    '''
    # Get gravity force in the global frame
    gravityForce = self.earthModel.getGravityForce(inertia, state)
    
    try:
        # Convert to local frame if in 6DoF simulation
        gravityForce = state.orientation.conjugate().rotate(gravityForce) # rotate into local frame when in 6DoF mode
    except AttributeError:
        pass # Don't do anything in 3DoF mode (No local frame exists)

    return ForceMomentSystem(gravityForce, inertia.CG)

def getWind(self, AGLElevation, time=None)

Returns wind in the North-East-Up (Y-X-Z) frame

Expand source code

def getWind(self, AGLElevation, time=None):
    ''' Returns wind in the North-East-Up (Y-X-Z) frame '''
    meanWindVel = self.meanWindModel.getMeanWind(AGLElevation)
    
    if self.turbulenceModel != None:
        turbWindVel = self.turbulenceModel.getTurbVelocity(AGLElevation, meanWindVel, time)
    else:
        turbWindVel = Vector(0,0,0)

    return meanWindVel, turbWindVel

class EnvironmentalConditions (ASLAltitude, Temp, Pressure, Density, DynamicViscosity, Wind, MeanWind, TurbWind)

EnvironmentalConditions(ASLAltitude, Temp, Pressure, Density, DynamicViscosity, Wind, MeanWind, TurbWind)

Ancestors

builtins.tuple

Instance variables

var ASLAltitude: Alias for field number 0
var Density: Alias for field number 3
var DynamicViscosity: Alias for field number 4
var MeanWind: Alias for field number 6
var Pressure: Alias for field number 2
var Temp: Alias for field number 1
var TurbWind: Alias for field number 7
var Wind: Alias for field number 5

class FlatEarth

Interface for all earth models

Expand source code

class FlatEarth(NoEarth):
    # GM = 398600.5 # km^3/s^2 Geocentric gravitational constant (Mass of atmosphere included)
    GM = 398600.5e9 # m^3/s^2

    def getGravityForce(self, inertia, state):
        ''' Returns a gravity force vector in the global frame '''
        gravityDirection = Vector(0, 0, -1) # Down in the launch tower frame

        distanceFromCenterOfEarth = state.position.Z + 6371000 # Adding earth's average radius
        gravityMagnitude = inertia.mass * self.GM / distanceFromCenterOfEarth**2 # Checked that this equation gives results similar to the 
            # US Standard Atmosphere - offset a bit depending on what earth radius is used. Current 6371 km gives slightly lower values that USTDA.

        return gravityDirection * gravityMagnitude

Ancestors

NoEarth
MAPLEAF.ENV.EarthModelling.EarthModel
abc.ABC

Class variables

var GM

Methods

def getGravityForce(self, inertia, state)

Returns a gravity force vector in the global frame

Expand source code

def getGravityForce(self, inertia, state):
    ''' Returns a gravity force vector in the global frame '''
    gravityDirection = Vector(0, 0, -1) # Down in the launch tower frame

    distanceFromCenterOfEarth = state.position.Z + 6371000 # Adding earth's average radius
    gravityMagnitude = inertia.mass * self.GM / distanceFromCenterOfEarth**2 # Checked that this equation gives results similar to the 
        # US Standard Atmosphere - offset a bit depending on what earth radius is used. Current 6371 km gives slightly lower values that USTDA.

    return gravityDirection * gravityMagnitude

Inherited members

NoEarth:

class Hellman (groundMeanWind=<MAPLEAF.Motion.CythonVector.Vector object>, altitudeLimit=1500, hellmanAlphaCoeff=0.14)

Uses the Hellman law to scale a ground velocity using a power law as the atmospheric boundary layer is exited. Low-altitude only. https://en.wikipedia.org/wiki/Wind_gradient -> Engineering Section.

Arguments

groundMeanWind: Wind Vector at ground level (~10m above surface), in m/s altitudeLimit: Altitude above which the power law transitions into a constant profile (m, AGL) hellmanAlphaCoeff: Alpha in the Hellman model (1/7 is a common first guess, but can vary widely depending on terrain)

Expand source code

class Hellman(MeanWindModel):
    '''
        Uses the Hellman law to scale a ground velocity using a power law as the atmospheric boundary layer is exited.
        Low-altitude only. https://en.wikipedia.org/wiki/Wind_gradient -> Engineering Section.

    '''
    def __init__(self, groundMeanWind=Vector(0,0,0), altitudeLimit=1500, hellmanAlphaCoeff=0.14):
        '''
            Arguments:
                groundMeanWind: Wind Vector at ground level (~10m above surface), in m/s
                altitudeLimit: Altitude above which the power law transitions into a constant profile (m, AGL)
                hellmanAlphaCoeff: Alpha in the Hellman model (1/7 is a common first guess, but can vary widely depending on terrain)
        '''
        self.groundMeanWind = groundMeanWind
        self.altitudeLimit = altitudeLimit
        self.hellmanAlphaCoeff = hellmanAlphaCoeff

    def getMeanWind(self, AGLAltitude):
        # Limit velocity scaling to up to the specified altitude limit
        HellmanAltitude = max(min(self.altitudeLimit, AGLAltitude), 10)
        # Assume initial winds come from a height of 10m
        return self.groundMeanWind * (HellmanAltitude/10)**self.hellmanAlphaCoeff

Ancestors

MAPLEAF.ENV.MeanWindModelling.MeanWindModel
abc.ABC

Methods

def getMeanWind(self, AGLAltitude)

Expand source code

def getMeanWind(self, AGLAltitude):
    # Limit velocity scaling to up to the specified altitude limit
    HellmanAltitude = max(min(self.altitudeLimit, AGLAltitude), 10)
    # Assume initial winds come from a height of 10m
    return self.groundMeanWind * (HellmanAltitude/10)**self.hellmanAlphaCoeff

class InterpolatedWind (windAltitudes=[], winds=[], windFilePath=None)

Defines the interface for Mean Wind Models

Arguments

windAltitudes: list of altitudes at which wind vectors will be provided (m AGL) winds: list of wind vectors, matching the ordering of windAltitudes (m/s)

windFilePath: string path to a text file containing a table, formatted similarly to the following: AGlAltitude(m) WindX(m/s) WindY(m/s) WindZ(m/s) a1 wx1 wy1 wz1 a2 wx2 wy2 wz2

Notes

Provide either windAltitudes AND winds, OR windFilePath. If all are provided, windAltitudes/winds will be used

Expand source code

class InterpolatedWind(MeanWindModel):
    def __init__(self, windAltitudes=[], winds=[], windFilePath=None):
        '''
            Arguments:
                windAltitudes: list of altitudes at which wind vectors will be provided (m AGL)
                winds: list of wind vectors, matching the ordering of windAltitudes (m/s)

                windFilePath: string path to a text file containing a table, formatted similarly to the following:
                    AGlAltitude(m) WindX(m/s) WindY(m/s) WindZ(m/s)
                    a1              wx1         wy1     wz1
                    a2              wx2         wy2     wz2

            Notes:
                Provide either windAltitudes AND winds, OR windFilePath. If all are provided, windAltitudes/winds will be used
        '''
        if len(windAltitudes) + len(winds) > 0:
            self.windAltitudes = windAltitudes
            self.winds = winds

        elif windFilePath != None:
            self._readCustomWindProfile(windFilePath)

        else:
            raise ValueError("Incomplete initialization info provided. Provide either windAltitudes AND winds, OR windFilePath")

    def _readCustomWindProfile(self, filePath):
        windProfile = np.loadtxt(filePath, skiprows=1)
        self.windAltitudes = windProfile[:, 0]
        self.winds = windProfile[:, 1:]

    def getMeanWind(self, AGLAltitude):
        return Vector(*linInterp(self.windAltitudes, self.winds, AGLAltitude))

Ancestors

MAPLEAF.ENV.MeanWindModelling.MeanWindModel
abc.ABC

Methods

def getMeanWind(self, AGLAltitude)

Expand source code

def getMeanWind(self, AGLAltitude):
    return Vector(*linInterp(self.windAltitudes, self.winds, AGLAltitude))

class LaunchRail (initialPosition, initialDirection, length, earthRotationRate=0)

Provide initial position and direction in global inertial frame

Expand source code

class LaunchRail():
    def __init__(self, initialPosition, initialDirection, length, earthRotationRate=0):
        ''' Provide initial position and direction in global inertial frame '''
        self.initialPosition = initialPosition # m
        self.initialDirection = initialDirection.normalize() # Ensure the direction is a unit vector
        self.length = length # m
        self.earthAngVel = earthRotationRate # rad/s, about z-axis

    def getPosition(self, time):
        ''' Launch Rail moves with the earth's surface '''
        rotationSoFar = Quaternion(
            axisOfRotation=Vector(0,0,1), 
            angle=(self.earthAngVel * time)
        )
        return rotationSoFar.rotate(self.initialPosition)

    def getDirection(self, time):
        ''' Launch Rail moves with the earth's surface '''
        rotationSoFar = Quaternion(
            axisOfRotation=Vector(0,0,1), 
            angle=(self.earthAngVel * time)
        )
        return rotationSoFar.rotate(self.initialDirection)

    def getVelocity(self, time):
        pos = self.getPosition(time)
        angVel = AngularVelocity(0, 0, self.earthAngVel)
        return angVel.crossProduct(pos)

    def applyLaunchTowerForce(self, state, time, unadjustedForce):
        '''
            If on launch tower, projects forces experienced onto the launch tower directions and sets Moments = 0
            Returns two values: onTower(Bool), adjustedForce(Motion.ForceMomentSystem)
        '''
        distanceTravelled = (state.position - self.getPosition(time)).length()
        if distanceTravelled < self.length:
            # Vehicle still on launch rail

            # Project total force onto the launch rail direction (dot product)
            adjustedForceMagnitude = unadjustedForce.force * self.getDirection(time)
            adjustedForceVector = self.initialDirection * adjustedForceMagnitude

            # No resultant moments while on the launch rail
            adjustedForce = ForceMomentSystem(adjustedForceVector)

            return adjustedForce
        else:
            # Vehicle has left the rail
            return unadjustedForce

    def applyLaunchTowerMotionConstraint(self, state, time):
        ''' If on launch tower, stops rocket from sliding off the bottom before engine is lit '''
        currPosition = self.getPosition(time)
        currDirection = self.getDirection(time)

        distanceTravelled = (state.position - currPosition) * currDirection
        if distanceTravelled < self.length:
            # Vehicle still on launch rail

            if distanceTravelled <= 0:
                # If vehicle has slid of the bottom of the rail
                # Reset its position to the bottom of the rail
                state.position = currPosition

                # Velocity adjustment
                velocityAlongRail = state.velocity * currDirection
                if velocityAlongRail < 0:
                    # If velocity along rail is also negative, set it to zero relative to the rail
                    state.velocity = self.getVelocity(time)

            return True, state
        else:
            # Vehicle has left the rail
            return False, state

Methods

def applyLaunchTowerForce(self, state, time, unadjustedForce)

If on launch tower, projects forces experienced onto the launch tower directions and sets Moments = 0 Returns two values: onTower(Bool), adjustedForce(Motion.ForceMomentSystem)

Expand source code

def applyLaunchTowerForce(self, state, time, unadjustedForce):
    '''
        If on launch tower, projects forces experienced onto the launch tower directions and sets Moments = 0
        Returns two values: onTower(Bool), adjustedForce(Motion.ForceMomentSystem)
    '''
    distanceTravelled = (state.position - self.getPosition(time)).length()
    if distanceTravelled < self.length:
        # Vehicle still on launch rail

        # Project total force onto the launch rail direction (dot product)
        adjustedForceMagnitude = unadjustedForce.force * self.getDirection(time)
        adjustedForceVector = self.initialDirection * adjustedForceMagnitude

        # No resultant moments while on the launch rail
        adjustedForce = ForceMomentSystem(adjustedForceVector)

        return adjustedForce
    else:
        # Vehicle has left the rail
        return unadjustedForce

def applyLaunchTowerMotionConstraint(self, state, time)

If on launch tower, stops rocket from sliding off the bottom before engine is lit

Expand source code

def applyLaunchTowerMotionConstraint(self, state, time):
    ''' If on launch tower, stops rocket from sliding off the bottom before engine is lit '''
    currPosition = self.getPosition(time)
    currDirection = self.getDirection(time)

    distanceTravelled = (state.position - currPosition) * currDirection
    if distanceTravelled < self.length:
        # Vehicle still on launch rail

        if distanceTravelled <= 0:
            # If vehicle has slid of the bottom of the rail
            # Reset its position to the bottom of the rail
            state.position = currPosition

            # Velocity adjustment
            velocityAlongRail = state.velocity * currDirection
            if velocityAlongRail < 0:
                # If velocity along rail is also negative, set it to zero relative to the rail
                state.velocity = self.getVelocity(time)

        return True, state
    else:
        # Vehicle has left the rail
        return False, state

def getDirection(self, time)

Launch Rail moves with the earth's surface

Expand source code

def getDirection(self, time):
    ''' Launch Rail moves with the earth's surface '''
    rotationSoFar = Quaternion(
        axisOfRotation=Vector(0,0,1), 
        angle=(self.earthAngVel * time)
    )
    return rotationSoFar.rotate(self.initialDirection)

def getPosition(self, time)

Launch Rail moves with the earth's surface

Expand source code

def getPosition(self, time):
    ''' Launch Rail moves with the earth's surface '''
    rotationSoFar = Quaternion(
        axisOfRotation=Vector(0,0,1), 
        angle=(self.earthAngVel * time)
    )
    return rotationSoFar.rotate(self.initialPosition)

def getVelocity(self, time)

Expand source code

def getVelocity(self, time):
    pos = self.getPosition(time)
    angVel = AngularVelocity(0, 0, self.earthAngVel)
    return angVel.crossProduct(pos)

class NoEarth

Interface for all earth models

Expand source code

class NoEarth(EarthModel):
    rotationRate = 0 # rad/s
    noRotation = Quaternion(1,0,0,0)
    noForce = Vector(0,0,0)

    def getGravityForce(self, inertia, state):
        return self.noForce

    def getInertialToENUFrameRotation(self, x, y, z):
        return self.noRotation

    def getAltitude(self, x, y, z):
        return z

    def convertIntoGlobalFrame(self, launchTowerFrameState, _, __):
        ''' Launch tower frame is global inertial frame '''
        return launchTowerFrameState

Ancestors

MAPLEAF.ENV.EarthModelling.EarthModel
abc.ABC

Subclasses

FlatEarth

Class variables

var noForce
var noRotation
var rotationRate

Methods

def convertIntoGlobalFrame(self, launchTowerFrameState, _, __)

Launch tower frame is global inertial frame

Expand source code

def convertIntoGlobalFrame(self, launchTowerFrameState, _, __):
    ''' Launch tower frame is global inertial frame '''
    return launchTowerFrameState

def getAltitude(self, x, y, z)

Given the aircraft coordinates in the global inertial frame, should return the altitude above SEA level (ASL) - to be used in wind calculations

Expand source code

def getAltitude(self, x, y, z):
    return z

def getGravityForce(self, inertia, state)

Return gravity force vector in the global inertial frame

Expand source code

def getGravityForce(self, inertia, state):
    return self.noForce

def getInertialToENUFrameRotation(self, x, y, z)

Return Quaternion defining the rotation from the global inertial frame to the North-East-Up frame defined on the surface under the air vehicle's position

Expand source code

def getInertialToENUFrameRotation(self, x, y, z):
    return self.noRotation

class PinkNoiseGenerator (alpha=1.6666666666666667, seed=None, nPoles=2, simulatedSamplingFreqHz=20)

Inputs

alpha - pink noise generated has a power spectrum which has a log-log slope of -alpha.
default value of 5/3 matches power spectrum slope of turbulent integral scales
seed - pink noise generated based on applying a filter to gaussian white noise specifying
the seed for the input random noise allows for repeatable sequences of values
nPoles - more poles include more low frequency content in the spectrum
parameter could be interpreted as similar to the turbulence length scale
can be varied in more sophisticated models
simulatedSamplingFreqHz - controls "sampling frequency" that the pink noise is interpreted as
also modulated the effective length scale of turbulence, together with nPoles

Methods

.getValue(time)
If the value is requested at time values that fall between sampling intervals, linearly interpolated values (from adjacent samples) are returned. If no time is provided, the sampling frequency is ignored and a new value is computed

Expand source code

class PinkNoiseGenerator():
    '''
        Inputs:
            
            * alpha - pink noise generated has a power spectrum which has a log-log slope of -alpha.  
                default value of 5/3 matches power spectrum slope of turbulent integral scales  
            * seed - pink noise generated based on applying a filter to gaussian white noise specifying  
                the seed for the input random noise allows for repeatable sequences of values  
            * nPoles - more poles include more low frequency content in the spectrum  
                parameter could be interpreted as similar to the turbulence length scale  
                can be varied in more sophisticated models  
            * simulatedSamplingFreqHz - controls "sampling frequency" that the pink noise is interpreted as  
                also modulated the effective length scale of turbulence, together with nPoles  
        Methods:
            
            * .getValue(time)  
                If the value is requested at time values that fall between sampling intervals, linearly 
                interpolated values (from adjacent samples) are returned.
                If no time is provided, the sampling frequency is ignored and a new value is computed  
    '''

    def __init__(self, alpha=5/3, seed=None, nPoles=2, simulatedSamplingFreqHz=20):
        self.alpha = alpha
        if seed != None:
            self.seed = seed
        else:
            self.seed = random.randrange(1000000)
        self.rng = random.Random(self.seed)
        # self.rng.seed(self.seed)
        self.nPoles = nPoles

        # Calculate coefficient values - openrocket documentation pg. 59/60 (Section 4.1.2 Wind Modeling)
        # Based on method from Kasdin (1995) "Discrete Simulation of Colored Noise and Stochastic Processes and 1/f^{alpha} Power Law Noise Generation"
        self.coeffs = [ 1, ]
        for i in range(nPoles):
            self.coeffs.append((i - alpha/2) * self.coeffs[i]/(i+1))
        self.coeffs.pop(0)

        # Overwriting will wrap around the array to avoid constantly shifting values
        self.nextIndexToOverwrite = 0

        # Fixed sampling frequency allows for achieving "mesh convergence" as the timestep is decreased
        self.lastValTime = 0
        self.timeStep = 1 / simulatedSamplingFreqHz # seconds

        # Fill lastValues
        self.lastValues = [ self.rng.gauss(0, 1) for i in range(nPoles) ]
        for i in range(10*nPoles):
            self.getValue()

    def getValue(self, time=None):
        '''
            Calculate coefficient values - openrocket documentation pg. 59/60 (Section 4.1.2 Wind Modeling)
            Based on method from Kasdin (1995) "Discrete Simulation of Colored Noise and Stochastic Processes and 1/f^{alpha} Power Law Noise Generation"
        '''
        # Compute any required new values
        while time == None or time > self.lastValTime:
            whiteNoise = self.rng.gauss(0, 1)

            for i in range(self.nPoles):
                whiteNoise -= self.coeffs[i]*self.lastValues[(self.nextIndexToOverwrite - i - 1) % self.nPoles]
            pinkNoiseSample = whiteNoise

            self.lastValues[self.nextIndexToOverwrite] = pinkNoiseSample
            self.nextIndexToOverwrite = (self.nextIndexToOverwrite + 1) % self.nPoles

            if time == None:
                return pinkNoiseSample

            self.lastValTime += self.timeStep
            if self.lastValTime == time:
                return pinkNoiseSample

        # Linearly interpolate value at desired time
        lastVal = self.lastValues[self.nextIndexToOverwrite - 1]
        secondLastVal = self.lastValues[self.nextIndexToOverwrite - 2]
        secondLastValTime = self.lastValTime - self.timeStep
        
        dValdt = (lastVal - secondLastVal) / self.timeStep
        return dValdt*(time-secondLastValTime) + secondLastVal

Methods

def getValue(self, time=None)

Calculate coefficient values - openrocket documentation pg. 59/60 (Section 4.1.2 Wind Modeling) Based on method from Kasdin (1995) "Discrete Simulation of Colored Noise and Stochastic Processes and 1/f^{alpha} Power Law Noise Generation"

Expand source code

def getValue(self, time=None):
    '''
        Calculate coefficient values - openrocket documentation pg. 59/60 (Section 4.1.2 Wind Modeling)
        Based on method from Kasdin (1995) "Discrete Simulation of Colored Noise and Stochastic Processes and 1/f^{alpha} Power Law Noise Generation"
    '''
    # Compute any required new values
    while time == None or time > self.lastValTime:
        whiteNoise = self.rng.gauss(0, 1)

        for i in range(self.nPoles):
            whiteNoise -= self.coeffs[i]*self.lastValues[(self.nextIndexToOverwrite - i - 1) % self.nPoles]
        pinkNoiseSample = whiteNoise

        self.lastValues[self.nextIndexToOverwrite] = pinkNoiseSample
        self.nextIndexToOverwrite = (self.nextIndexToOverwrite + 1) % self.nPoles

        if time == None:
            return pinkNoiseSample

        self.lastValTime += self.timeStep
        if self.lastValTime == time:
            return pinkNoiseSample

    # Linearly interpolate value at desired time
    lastVal = self.lastValues[self.nextIndexToOverwrite - 1]
    secondLastVal = self.lastValues[self.nextIndexToOverwrite - 2]
    secondLastValTime = self.lastValTime - self.timeStep
    
    dValdt = (lastVal - secondLastVal) / self.timeStep
    return dValdt*(time-secondLastValTime) + secondLastVal

class SphericalEarth

Models a non-rotating, uniform, spherical earth

Expand source code

class SphericalEarth(EarthModel):
    ''' Models a non-rotating, uniform, spherical earth '''
    radius = 6371007.1809 # m - from NESC-RP-12-007770, Volume II
    rotationRate = 7.292115e-5 # rad/s angular velocity - from WGS84 model. Defined WRT stars (not our sun)
    # GM = 398600.5 # km^3/s^2 Geocentric gravitational constant (Mass of atmosphere included)
    GM = 3.98600436e14 # m^3/s^2- from NESC-RP-12-007770, Volume II

    def geodeticToCartesian(self, lat, lon, height, timeOffset=0):
        '''
            Converts from geodetic coordinates to cartesian coordinates, assuming the surface is a sphere.
            
            Resulting Z axis goes through north pole, X axis goes through prime meridian at the equator, and 
                the Y axis is perpendicular to both in such a way as to form a right-handed coordinate system

            Inputs:
                lat:    (float) latitude - degrees
                lon:    (float) longitude - degrees
                height: (float) ASL (usually m)

            Returns:
                x:      (float) meters, relative to center of sphere
                y:      ' '
                z:      ' '
        '''
        lat = radians(lat)
        lon = radians(lon)

        # Adjust for rotation rate
        lon += self.rotationRate * timeOffset

        radius = self.radius + height

        x = radius * cos(lat) * cos(lon)
        y = radius * cos(lat) * sin(lon)
        z = radius * sin(lat)
        return x, y, z

    def cartesianToGeodetic(self, x, y, z, timeOffset=0):
        ''' Convert from cartesian to lat/lon/height coordinates, where height is ASL '''
        p = sqrt(x*x + y*y)
        
        lon = atan2(y, x)

        # Adjust for rotation rate
        lon -= self.rotationRate * timeOffset

        lat = atan2(z, p)
        height = sqrt(x*x + y*y + z*z) - self.radius
        
        return degrees(lat), degrees(lon), height

    def getGravityForce(self, inertia, state):
        ''' Returns a gravity force vector in the global frame '''
        gravityDirection = -(state.position.normalize()) # Gravity pulls toward the earth's center of mass, which is the origin
        
        distanceFromCenterOfEarth = state.position.length()
        gravityMagnitude = inertia.mass * self.GM / distanceFromCenterOfEarth**2 # Checked that this equation gives results similar to the 
            # US Standard Atmosphere - offset a bit depending on what earth radius is used. Current 6371 km gives slightly lower values that USTDA.

        return gravityDirection * gravityMagnitude

    def getInertialToENUFrameRotation(self, x, y, z):
        ''' Returns a Quaternion that defines the rotation b/w the global inertia frame and the local, surface-normal North-East-Up (y-x-z) frame '''
        # Time offset not necessary, since the x, y, z coordinates are inertial, and we are not 
            # interested in finding a particular location on the surface of the earth,
            # we are just interested in finding the surface normal under a x/y/z location. 
            # The result is independent of the earth's rotation.
        lat, lon, height = self.cartesianToGeodetic(x, y, z)

        # Rotation between inertial frame and ENU frame will be composed in two steps
        # Step 1: Rotate about inertial/initial Z-axis to Y-axis in such a way, that after rotation 2, 
            # it will be pointing North
        rot1Angle = radians(lon + 90)
        rot1 = Quaternion(axisOfRotation=Vector(0,0,1), angle=rot1Angle)
        
        # Step 2: Rotate about the local x-axis to match the latitude
        rot2Angle = radians(90 - lat)
        rot2 = Quaternion(axisOfRotation=Vector(1,0,0), angle=rot2Angle)

        # Sequentially combine rotations by multiplying them
        return rot1 * rot2

    def getAltitude(self, x, y, z):
        # Height is independent of the earth's rotation
        _, __, height = self.cartesianToGeodetic(x, y, z)
        return height

    def convertIntoGlobalFrame(self, launchTowerFrameState, lat, lon):
        ### Position ###
        height = launchTowerFrameState.position.Z # ASL altitude
        globalFramePosition = Vector(*self.geodeticToCartesian(lat, lon, height))
        
        ### Velocity ###
        # Find orientation of launch tower frame relative to global frame
        launchTowerToGlobalFrame = self.getInertialToENUFrameRotation(*globalFramePosition)

        # Rotate velocity accordingly
        rotatedVelocity = launchTowerToGlobalFrame.rotate(launchTowerFrameState.velocity)
        # Add earth's surface velocity
        earthAngVel = Vector(0, 0, self.rotationRate)
        velocityDueToEarthRotation = earthAngVel.crossProduct(globalFramePosition)
        globalFrameVelocity = rotatedVelocity + velocityDueToEarthRotation

        try:  # 6DoF Case           
            ### Orientation ###
            globalFrameOrientation = launchTowerToGlobalFrame * launchTowerFrameState.orientation

            ### Angular Velocity ###
            # Angular velocity is defined in the vehicle's local frame, so the conversion needs to go the other way
            earthAngVel_RocketFrame = globalFrameOrientation.conjugate().rotate(earthAngVel)
            localFrame_adjustedAngVel = launchTowerFrameState.angularVelocity + earthAngVel_RocketFrame

            return RigidBodyState( 
                globalFramePosition, 
                globalFrameVelocity, 
                globalFrameOrientation, 
                localFrame_adjustedAngVel
                )
        
        except AttributeError: # 3DoF Case
            return RigidBodyState_3DoF(globalFramePosition, globalFrameVelocity)

Ancestors

MAPLEAF.ENV.EarthModelling.EarthModel
abc.ABC

Subclasses

WGS84

Class variables

var GM
var radius
var rotationRate

Methods

def cartesianToGeodetic(self, x, y, z, timeOffset=0)

Convert from cartesian to lat/lon/height coordinates, where height is ASL

Expand source code

def cartesianToGeodetic(self, x, y, z, timeOffset=0):
    ''' Convert from cartesian to lat/lon/height coordinates, where height is ASL '''
    p = sqrt(x*x + y*y)
    
    lon = atan2(y, x)

    # Adjust for rotation rate
    lon -= self.rotationRate * timeOffset

    lat = atan2(z, p)
    height = sqrt(x*x + y*y + z*z) - self.radius
    
    return degrees(lat), degrees(lon), height

def convertIntoGlobalFrame(self, launchTowerFrameState, lat, lon)

Should take a RigidBodyState defined in the launch tower frame (fixed to the earth's surface), where position.Z has been redefined relative to sea level (instead of ground-level), and convert it into the global inertial frame. Exception is the Angular Velocity, since it is defined in the vehicle's local frame. If the earth model is rotating, add the earth's angular velocity to it (after conversion into the local frame)

Expand source code

def convertIntoGlobalFrame(self, launchTowerFrameState, lat, lon):
    ### Position ###
    height = launchTowerFrameState.position.Z # ASL altitude
    globalFramePosition = Vector(*self.geodeticToCartesian(lat, lon, height))
    
    ### Velocity ###
    # Find orientation of launch tower frame relative to global frame
    launchTowerToGlobalFrame = self.getInertialToENUFrameRotation(*globalFramePosition)

    # Rotate velocity accordingly
    rotatedVelocity = launchTowerToGlobalFrame.rotate(launchTowerFrameState.velocity)
    # Add earth's surface velocity
    earthAngVel = Vector(0, 0, self.rotationRate)
    velocityDueToEarthRotation = earthAngVel.crossProduct(globalFramePosition)
    globalFrameVelocity = rotatedVelocity + velocityDueToEarthRotation

    try:  # 6DoF Case           
        ### Orientation ###
        globalFrameOrientation = launchTowerToGlobalFrame * launchTowerFrameState.orientation

        ### Angular Velocity ###
        # Angular velocity is defined in the vehicle's local frame, so the conversion needs to go the other way
        earthAngVel_RocketFrame = globalFrameOrientation.conjugate().rotate(earthAngVel)
        localFrame_adjustedAngVel = launchTowerFrameState.angularVelocity + earthAngVel_RocketFrame

        return RigidBodyState( 
            globalFramePosition, 
            globalFrameVelocity, 
            globalFrameOrientation, 
            localFrame_adjustedAngVel
            )
    
    except AttributeError: # 3DoF Case
        return RigidBodyState_3DoF(globalFramePosition, globalFrameVelocity)

def geodeticToCartesian(self, lat, lon, height, timeOffset=0)

Converts from geodetic coordinates to cartesian coordinates, assuming the surface is a sphere.

Resulting Z axis goes through north pole, X axis goes through prime meridian at the equator, and the Y axis is perpendicular to both in such a way as to form a right-handed coordinate system

Inputs

lat: (float) latitude - degrees lon: (float) longitude - degrees height: (float) ASL (usually m)

Returns

x: (float) meters, relative to center of sphere
y: ' '
z: ' '

Expand source code

def geodeticToCartesian(self, lat, lon, height, timeOffset=0):
    '''
        Converts from geodetic coordinates to cartesian coordinates, assuming the surface is a sphere.
        
        Resulting Z axis goes through north pole, X axis goes through prime meridian at the equator, and 
            the Y axis is perpendicular to both in such a way as to form a right-handed coordinate system

        Inputs:
            lat:    (float) latitude - degrees
            lon:    (float) longitude - degrees
            height: (float) ASL (usually m)

        Returns:
            x:      (float) meters, relative to center of sphere
            y:      ' '
            z:      ' '
    '''
    lat = radians(lat)
    lon = radians(lon)

    # Adjust for rotation rate
    lon += self.rotationRate * timeOffset

    radius = self.radius + height

    x = radius * cos(lat) * cos(lon)
    y = radius * cos(lat) * sin(lon)
    z = radius * sin(lat)
    return x, y, z

def getAltitude(self, x, y, z)

Given the aircraft coordinates in the global inertial frame, should return the altitude above SEA level (ASL) - to be used in wind calculations

Expand source code

def getAltitude(self, x, y, z):
    # Height is independent of the earth's rotation
    _, __, height = self.cartesianToGeodetic(x, y, z)
    return height

def getGravityForce(self, inertia, state)

Returns a gravity force vector in the global frame

Expand source code

def getGravityForce(self, inertia, state):
    ''' Returns a gravity force vector in the global frame '''
    gravityDirection = -(state.position.normalize()) # Gravity pulls toward the earth's center of mass, which is the origin
    
    distanceFromCenterOfEarth = state.position.length()
    gravityMagnitude = inertia.mass * self.GM / distanceFromCenterOfEarth**2 # Checked that this equation gives results similar to the 
        # US Standard Atmosphere - offset a bit depending on what earth radius is used. Current 6371 km gives slightly lower values that USTDA.

    return gravityDirection * gravityMagnitude

def getInertialToENUFrameRotation(self, x, y, z)

Returns a Quaternion that defines the rotation b/w the global inertia frame and the local, surface-normal North-East-Up (y-x-z) frame

Expand source code

def getInertialToENUFrameRotation(self, x, y, z):
    ''' Returns a Quaternion that defines the rotation b/w the global inertia frame and the local, surface-normal North-East-Up (y-x-z) frame '''
    # Time offset not necessary, since the x, y, z coordinates are inertial, and we are not 
        # interested in finding a particular location on the surface of the earth,
        # we are just interested in finding the surface normal under a x/y/z location. 
        # The result is independent of the earth's rotation.
    lat, lon, height = self.cartesianToGeodetic(x, y, z)

    # Rotation between inertial frame and ENU frame will be composed in two steps
    # Step 1: Rotate about inertial/initial Z-axis to Y-axis in such a way, that after rotation 2, 
        # it will be pointing North
    rot1Angle = radians(lon + 90)
    rot1 = Quaternion(axisOfRotation=Vector(0,0,1), angle=rot1Angle)
    
    # Step 2: Rotate about the local x-axis to match the latitude
    rot2Angle = radians(90 - lat)
    rot2 = Quaternion(axisOfRotation=Vector(1,0,0), angle=rot2Angle)

    # Sequentially combine rotations by multiplying them
    return rot1 * rot2

class TabulatedAtmosphere (filePath)

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)

Expand source code

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)

Ancestors

MAPLEAF.ENV.AtmosphereModelling.AtmosphericModel
abc.ABC

Methods

def getAirProperties(self, ASLElevation)

Return an iterable containing: temp(K), static pressure (Pa), density (kg/m^3), dynamic viscosity (Pa*s), in that order

Expand source code

def getAirProperties(self, ASLElevation, _=None):
    return linInterp(self.keys, self.values, ASLElevation)

class USStandardAtmosphere

Provides atmospheric properties calculated according to the model here: https://nebula.wsimg.com/ab321c1edd4fa69eaa94b5e8e769b113?AccessKeyId=AF1D67CEBF3A194F66A3&disposition=0&alloworigin=1

Expand source code

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 ]

Ancestors

MAPLEAF.ENV.AtmosphereModelling.AtmosphericModel
abc.ABC

Class variables

var G
var M
var R
var T0
var baseHeights
var basePressures
var baseTemps
var dt_dh
var earthRadius
var p0

Methods

def getAirProperties(self, ASLElevation, time)

Return an iterable containing: temp(K), static pressure (Pa), density (kg/m^3), dynamic viscosity (Pa*s), in that order

Expand source code

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 ]

class WGS84

Models a rotating, ellipsoid earth, with non-uniform gravity Inherits the getInertialToENUFrameRotation function from SphericalEarth, otherwise overrides everything

Expand source code

class WGS84(SphericalEarth):
    ''' 
        Models a rotating, ellipsoid earth, with non-uniform gravity 
        Inherits the getInertialToENUFrameRotation function from SphericalEarth, otherwise overrides everything
    '''
    # Set up/calculate the WGS84 Ellipsoid parameters and derived values
    a = 6378137             # m, semi-major axis
    f = 1/298.257223563     # Flattening (derived from C_20)
    b = a * (1-f)           # m, semi-minor axis = 6 356 752.314 140
    e2 = f*(2-f)            # First eccentricity squared (e^2) = 6.694 379 990 14 e-3
    eP2 = e2 / ((1-f)**2)   # Second eccentricity squared (e'^2) = 6.739 496 742 28 e-3s
    eP = sqrt(eP2)          # Second eccentricity (e')

    def __init__(self):
        # Read gravity coefficients from table
        # Columns are: Degree Order C S (Where C and S are the coefficients)
        coeffPath = getAbsoluteFilePath('./MAPLEAF/ENV/sphericalHarmonicGravityCoeffs.txt')
        gravityCoeffs = np.loadtxt(coeffPath, skiprows=2)
        
        # Convert that table into nested dictionary form, to make C and S coefficients easily accessible by degree and order
        Ccoeffs = {}
        Scoeffs = {}
        for row in gravityCoeffs:
            degree, order, C, S, = row
            
            if degree not in Ccoeffs:
                # Create subdictionary for each degree if this is the first time we're encountering that degree
                Ccoeffs[degree] = {}
                Scoeffs[degree] = {}

            # Save coefficient from this row
            Ccoeffs[degree][order] = C
            Scoeffs[degree][order] = S

        # Save results
        # Now C_{2,3} is accessible as self.C[2][3], same for self.S
        self.C = Ccoeffs
        self.S = Scoeffs

    def geodeticToCartesian(self, lat, lon, height, timeOffset=0):
        '''        
            Converts from geodetic coordinates to cartesian coordinates, assuming the surface is a sphere.
            
            Resulting Z axis goes through north pole, X axis goes through prime meridian at the equator, and 
                the Y axis is perpendicular to both in such a way as to form a right-handed coordinate system
            
            Method adapted from: https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

            Note: In other references, it is common to set \phi = lat, \lamda = lon

            Inputs:
                lat:        (float) latitude - degrees
                lon:        (float) longitude - degrees
                height:     (float) ASL - m
                timeOffset  (float) Offset added to results based on earth's rotation rate.
                                At timeOffset=0, this inertial frame is aligned with the standard earth-centered-earth-fixed (ECEF) frame

            Returns:
                x:      (float) meters, relative to center of ellipsoid
                y:      ' '
                z:      ' '
        '''
        lat = radians(lat)
        lon = radians(lon)

        # Adjust for earth's rate of rotation
        lon += self.rotationRate * timeOffset

        N = self.a / sqrt(1 - self.e2*sin(lat)**2)

        x = (N + height)*cos(lat)*cos(lon)
        y = (N + height)*cos(lat)*sin(lon)
        z = ((self.b**2 / self.a**2)*N + height) * sin(lat)

        return x, y, z

    def cartesianToGeodetic(self, x, y, z, timeOffset=0):
        '''
            Inputs in m, sec
            Convert from earth-centered inertial coordinates to WGS84 geodetic coordinates
            Method adapted from: https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

            At timeOffset=0, the inertial frame is aligned with the standard earth-centered-earth-fixed (ECEF) frame

            Returns latitude, longitude, height in degrees, m
        '''
        lon = atan2(y, x) # Longitude is easy to find

        # Adjust longitude for earth's rotation
        lon -= self.rotationRate * timeOffset

        ### Compute latitude and height normal to surface w/ Newton-Raphson method ###
        p = sqrt(x*x + y*y)
        Kold = 1/(1 - self.e2)
        Korig = Kold

        while True:
            ci = ((p*p + (1-self.e2)*z*z * Kold**2)**(1.5)) / (self.a * self.e2)
            Knew = 1 + (p*p + (1 - self.e2)*z*z * Kold**3) / (ci - p*p)
            
            # End iterations once delta per iteration < 1e-10
            if abs(Knew - Kold) < 1e-10:
                break

            Kold = Knew

        k = Knew

        lat = atan2(k*z, p)

        h = (1/self.e2)*(1/k - 1/Korig)*sqrt(p*p + z*z*k*k)

        return degrees(lat), degrees(lon), h

    # Re-use conversions to/from spherical coordinates from SphericalEarth
        # Useful for Spherical Harmonic Gravity model calculations, which are performed in spherical coordinates
    cartesianToSpherical = SphericalEarth.cartesianToGeodetic
    sphericalToCartesian = SphericalEarth.geodeticToCartesian

    def getGravityForce(self, inertia, state):
        ''' 
            Get gravity using the J2 model (just a single spherical harmonic coefficient accounting for the earth's oblateness) 
            Method from NESC-RP-12-00770, Volume II, pg.51, Eqn 29-31
        '''
        # J2 Model just uses the first spherical harmonic coefficient
        J2 = 0.00108262982
        re = self.a # Equatorial radius (Semi-major axis of ellipsoid)

        x, y, z = state.position[0], state.position[1], state.position[2]
        r = state.position.length()
        
        mu = self.GM

        # Precompute common parts of equation
        frac = (3*J2 * re*re) / (2 * r**4)
        multiplier = -mu*inertia.mass / r**3

        xForce = x * multiplier * (1 - frac*(5*z*z - r*r))
        yForce = y * multiplier * (1 - frac*(5*z*z - r*r))
        zForce = z * multiplier * (1 - frac*(5*z*z - 3*r*r))

        return Vector(xForce, yForce, zForce)

Ancestors

SphericalEarth
MAPLEAF.ENV.EarthModelling.EarthModel
abc.ABC

Class variables

var a
var b
var e2
var eP
var eP2
var f

Methods

def cartesianToGeodetic(self, x, y, z, timeOffset=0)

Inputs in m, sec Convert from earth-centered inertial coordinates to WGS84 geodetic coordinates Method adapted from: https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

At timeOffset=0, the inertial frame is aligned with the standard earth-centered-earth-fixed (ECEF) frame

Returns latitude, longitude, height in degrees, m

Expand source code

def cartesianToGeodetic(self, x, y, z, timeOffset=0):
    '''
        Inputs in m, sec
        Convert from earth-centered inertial coordinates to WGS84 geodetic coordinates
        Method adapted from: https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

        At timeOffset=0, the inertial frame is aligned with the standard earth-centered-earth-fixed (ECEF) frame

        Returns latitude, longitude, height in degrees, m
    '''
    lon = atan2(y, x) # Longitude is easy to find

    # Adjust longitude for earth's rotation
    lon -= self.rotationRate * timeOffset

    ### Compute latitude and height normal to surface w/ Newton-Raphson method ###
    p = sqrt(x*x + y*y)
    Kold = 1/(1 - self.e2)
    Korig = Kold

    while True:
        ci = ((p*p + (1-self.e2)*z*z * Kold**2)**(1.5)) / (self.a * self.e2)
        Knew = 1 + (p*p + (1 - self.e2)*z*z * Kold**3) / (ci - p*p)
        
        # End iterations once delta per iteration < 1e-10
        if abs(Knew - Kold) < 1e-10:
            break

        Kold = Knew

    k = Knew

    lat = atan2(k*z, p)

    h = (1/self.e2)*(1/k - 1/Korig)*sqrt(p*p + z*z*k*k)

    return degrees(lat), degrees(lon), h

# Re-use conversions to/from spherical coordinates from SphericalEarth
    # Useful for Spherical Harmonic Gravity model calculations, which are performed in spherical coordinates

def cartesianToSpherical(self, x, y, z, timeOffset=0)

Convert from cartesian to lat/lon/height coordinates, where height is ASL

Expand source code

def cartesianToGeodetic(self, x, y, z, timeOffset=0):
    ''' Convert from cartesian to lat/lon/height coordinates, where height is ASL '''
    p = sqrt(x*x + y*y)
    
    lon = atan2(y, x)

    # Adjust for rotation rate
    lon -= self.rotationRate * timeOffset

    lat = atan2(z, p)
    height = sqrt(x*x + y*y + z*z) - self.radius
    
    return degrees(lat), degrees(lon), height

def geodeticToCartesian(self, lat, lon, height, timeOffset=0)

Converts from geodetic coordinates to cartesian coordinates, assuming the surface is a sphere.

Resulting Z axis goes through north pole, X axis goes through prime meridian at the equator, and the Y axis is perpendicular to both in such a way as to form a right-handed coordinate system

Method adapted from: https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

Note: In other references, it is common to set \phi = lat, \lamda = lon

Inputs

lat: (float) latitude - degrees lon: (float) longitude - degrees height: (float) ASL - m timeOffset (float) Offset added to results based on earth's rotation rate. At timeOffset=0, this inertial frame is aligned with the standard earth-centered-earth-fixed (ECEF) frame

Returns

x: (float) meters, relative to center of ellipsoid
y: ' '
z: ' '

Expand source code

def geodeticToCartesian(self, lat, lon, height, timeOffset=0):
    '''        
        Converts from geodetic coordinates to cartesian coordinates, assuming the surface is a sphere.
        
        Resulting Z axis goes through north pole, X axis goes through prime meridian at the equator, and 
            the Y axis is perpendicular to both in such a way as to form a right-handed coordinate system
        
        Method adapted from: https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

        Note: In other references, it is common to set \phi = lat, \lamda = lon

        Inputs:
            lat:        (float) latitude - degrees
            lon:        (float) longitude - degrees
            height:     (float) ASL - m
            timeOffset  (float) Offset added to results based on earth's rotation rate.
                            At timeOffset=0, this inertial frame is aligned with the standard earth-centered-earth-fixed (ECEF) frame

        Returns:
            x:      (float) meters, relative to center of ellipsoid
            y:      ' '
            z:      ' '
    '''
    lat = radians(lat)
    lon = radians(lon)

    # Adjust for earth's rate of rotation
    lon += self.rotationRate * timeOffset

    N = self.a / sqrt(1 - self.e2*sin(lat)**2)

    x = (N + height)*cos(lat)*cos(lon)
    y = (N + height)*cos(lat)*sin(lon)
    z = ((self.b**2 / self.a**2)*N + height) * sin(lat)

    return x, y, z

def getGravityForce(self, inertia, state)

Get gravity using the J2 model (just a single spherical harmonic coefficient accounting for the earth's oblateness) Method from NESC-RP-12-00770, Volume II, pg.51, Eqn 29-31

Expand source code

def getGravityForce(self, inertia, state):
    ''' 
        Get gravity using the J2 model (just a single spherical harmonic coefficient accounting for the earth's oblateness) 
        Method from NESC-RP-12-00770, Volume II, pg.51, Eqn 29-31
    '''
    # J2 Model just uses the first spherical harmonic coefficient
    J2 = 0.00108262982
    re = self.a # Equatorial radius (Semi-major axis of ellipsoid)

    x, y, z = state.position[0], state.position[1], state.position[2]
    r = state.position.length()
    
    mu = self.GM

    # Precompute common parts of equation
    frac = (3*J2 * re*re) / (2 * r**4)
    multiplier = -mu*inertia.mass / r**3

    xForce = x * multiplier * (1 - frac*(5*z*z - r*r))
    yForce = y * multiplier * (1 - frac*(5*z*z - r*r))
    zForce = z * multiplier * (1 - frac*(5*z*z - 3*r*r))

    return Vector(xForce, yForce, zForce)

def sphericalToCartesian(self, lat, lon, height, timeOffset=0)

Converts from geodetic coordinates to cartesian coordinates, assuming the surface is a sphere.

Resulting Z axis goes through north pole, X axis goes through prime meridian at the equator, and the Y axis is perpendicular to both in such a way as to form a right-handed coordinate system

Inputs

lat: (float) latitude - degrees lon: (float) longitude - degrees height: (float) ASL (usually m)

Returns

x: (float) meters, relative to center of sphere
y: ' '
z: ' '

Expand source code

def geodeticToCartesian(self, lat, lon, height, timeOffset=0):
    '''
        Converts from geodetic coordinates to cartesian coordinates, assuming the surface is a sphere.
        
        Resulting Z axis goes through north pole, X axis goes through prime meridian at the equator, and 
            the Y axis is perpendicular to both in such a way as to form a right-handed coordinate system

        Inputs:
            lat:    (float) latitude - degrees
            lon:    (float) longitude - degrees
            height: (float) ASL (usually m)

        Returns:
            x:      (float) meters, relative to center of sphere
            y:      ' '
            z:      ' '
    '''
    lat = radians(lat)
    lon = radians(lon)

    # Adjust for rotation rate
    lon += self.rotationRate * timeOffset

    radius = self.radius + height

    x = radius * cos(lat) * cos(lon)
    y = radius * cos(lat) * sin(lon)
    z = radius * sin(lat)
    return x, y, z

Inherited members

SphericalEarth: