Module `MAPLEAF.Rocket`

Code specific to modelling Rockets.
Main class is Rocket.

Other files define models for specific rocket components or stages.
In general the responsibility of rocket component classes is to model:

The forces/moments (aerodynamic or otherwise (excluding gravitational forces)) applied to the rocket by that rocket component
The inertia of that rocket component

The interface all rocket components are expected to implement is defined by the abstract base class RocketComponent. All rocket components inherit from this base class.

MAPLEAF.Rocket relies on (in order of coupling) MAPLEAF.Motion, MAPLEAF.ENV, MAPLEAF.IO, and MAPLEAF.GNC

Approximate Material Roughnesses

Material	Roughness ( $\mu m$ )
Mirror:	0
Glass:	0.1
Finished/Polished surface:	0.5
Aircraft-type sheet-metal surface:	2
Optimum paint-sprayed surfaces:	5
Planed wooden boards:	15
Paint in aircraft mass production:	20
Steel plating: bare:	50
Smooth cement:	50
Surface with asphalt-type coating:	100
Dip-galvanized metal surface:	150
Incorrectly sprayed paint:	200
Natural cast-iron surface:	250
Raw wooden boards:	500
Average concrete:	1000

(From Barrowman, 1967, Table 4-1)

Expand source code

r'''
Code specific to modelling Rockets.  
Main class is `Rocket`.

Other files define models for specific rocket components or stages.  
In general the responsibility of rocket component classes is to model:

1. The forces/moments (aerodynamic or otherwise (excluding gravitational forces)) applied to the rocket by that rocket component
2. The inertia of that rocket component

The interface all rocket components are expected to implement is defined by the abstract base class `MAPLEAF.Rocket.RocketComponent`. All rocket components inherit from this base class.

.. image:: https://airandspace.si.edu/sites/default/files/images/NASAJSC2002-01598h.jpg

MAPLEAF.Rocket relies on (in order of coupling) `MAPLEAF.Motion`, `MAPLEAF.ENV`, `MAPLEAF.IO`, and `MAPLEAF.GNC`

### Approximate Material Roughnesses
| Material                           | Roughness ( \(\mu m\) )                                        |
|------------------------------------|-------------------------------------------------------|
| Mirror:                            | 0                                                     |
| Glass:                             | 0.1                                                   |
| Finished/Polished surface:         | 0.5                                                   |
| Aircraft-type sheet-metal surface: | 2                                                     |
| Optimum paint-sprayed surfaces:    | 5                                                     |
| Planed wooden boards:              | 15                                                    |
| Paint in aircraft mass production: | 20                                                    |
| Steel plating: bare:               | 50                                                    |
| Smooth cement:                     | 50                                                    |
| Surface with asphalt-type coating: | 100                                                   |
| Dip-galvanized metal surface:      | 150                                                   |
| Incorrectly sprayed paint:         | 200                                                   |
| Natural cast-iron surface:         | 250                                                   |
| Raw wooden boards:                 | 500                                                   |
| Average concrete:                  | 1000                                                  |

(From Barrowman, 1967, Table 4-1)  

'''
# Make the classes in all submodules importable directly from MAPLEAF.Rocket
from .RocketComponents import *
from .simEventDetector import *

from .boatTail import *
from .bodyTube import *
from .Fins import *
from .Propulsion import *
from .sampleStatefulRocketComponent import *
from .Recovery import *
from .noseCone import *

from .RocketComponentFactory import * # Must be imported after the components above, the factory relies on them

from .stage import *
from .rocket import *

subModules = [ RocketComponents, simEventDetector, boatTail, bodyTube, Fins, Propulsion, sampleStatefulRocketComponent, Recovery, noseCone, stage, rocket, RocketComponentFactory ]

__all__ = [ ]

for subModule in subModules:
    __all__ += subModule.__all__

Sub-modules

MAPLEAF.Rocket.AeroFunctions: Functions to calculate parameters relevant to aerodynamic calculations - Mach/Reynolds numbers, AOA, local frame air velocity etc …
MAPLEAF.Rocket.CompositeObject: Both Stage and Rocket objects inherit from MAPLEAF.Rocket.CompositeObject. It implements functionality to add forces and inertias …
MAPLEAF.Rocket.CythonFinFunctions: Cython functions to speed up the FinSet aerodynamic model
MAPLEAF.Rocket.Fins
MAPLEAF.Rocket.Propulsion
MAPLEAF.Rocket.Recovery
MAPLEAF.Rocket.RocketComponentFactory: In charge of initializing rocket components. Add new components to stringNameToClassMap to make them instantiate themselves when included in a Rocket
MAPLEAF.Rocket.RocketComponents: Contains interface definitions (RocketComponent/BodyComponent) and Base classes for all RocketComponents (FixedMass),
as well as some simple …
MAPLEAF.Rocket.boatTail
MAPLEAF.Rocket.bodyTube
MAPLEAF.Rocket.noseCone
MAPLEAF.Rocket.rocket: Rocket ties together the code in MAPLEAF.GNC, MAPLEAF.Rocket, MAPLEAF.Motion, and MAPLEAF.ENV` to simulate the flight of a single rocket or …
MAPLEAF.Rocket.sampleStatefulRocketComponent
MAPLEAF.Rocket.simEventDetector: Generalized event detector (Apogee, motor burnout etc…) for a single vehicle.
Used to trigger stage separations and recovery system deployments.
MAPLEAF.Rocket.stage

Functions

def initializeForceLogging(component, subDictPath, rocket)

Expand source code

def initializeForceLogging(component, subDictPath, rocket):
    if rocket.derivativeEvaluationLog is not None:
        componentName = subDictPath[subDictPath.index(".")+1:]
        zeroVector = Vector(0,0,0)
        component.forcesLog = rocket.derivativeEvaluationLog.addColumn(componentName + "F(N)", zeroVector)
        component.momentsLog = rocket.derivativeEvaluationLog.addColumn(componentName + "M(Nm)", zeroVector)

def rocketComponentFactory(subDictPath, rocket, stage)

Initializes a rocket component based on the stringNameToClassMap

Inputs

subDictPath: (string) Path to subDict in simulation definition, like "Rocket.Stage1.Nosecone" rocket: (Rocket) that the component is a part of stage: (Stage) That the component is a part of Also uses the stringNameToClassMap dictionary

Expand source code

def rocketComponentFactory(subDictPath, rocket, stage):
    """
        Initializes a rocket component based on the stringNameToClassMap
        Inputs:
            subDictPath:        (string) Path to subDict in simulation definition, like "Rocket.Stage1.Nosecone"
            rocket:             (Rocket) that the component is a part of
            stage:              (Stage) That the component is a part of
        Also uses the stringNameToClassMap dictionary
    """       
    # Create SubDictReader for the rocket component's dictionary
    componentDictReader = SubDictReader(subDictPath, rocket.simDefinition)

    # Figure out which class to initialize
    className = componentDictReader.getString("class")
    referencedClass = stringNameToClassMap[className]
    
    # Initialize the rocket component
    newComponent = referencedClass(componentDictReader, rocket, stage)

    # Initialize logging component forces (if desired)
    initializeForceLogging(newComponent, subDictPath, rocket)

    return newComponent

Classes

class AeroDamping (componentDictReader, rocket, stage)

A zero-inertia component with constant aerodynamic damping coefficients

Expand source code

class AeroDamping(AeroForce):
    ''' A zero-inertia component with constant aerodynamic damping coefficients '''

    position = Vector(0,0,0)

    def __init__(self, componentDictReader, rocket, stage):
        self.componentDictReader = componentDictReader
        self.rocket = rocket
        self.stage = stage
        self.name = componentDictReader.getDictName()

        self.Aref = componentDictReader.getFloat("Aref")
        self.Lref = componentDictReader.getFloat("Lref")

        self.zDampingCoeffs = componentDictReader.getVector("zDampingCoeffs")
        self.yDampingCoeffs = componentDictReader.getVector("yDampingCoeffs")
        self.xDampingCoeffs = componentDictReader.getVector("xDampingCoeffs")
    
    def getAppliedForce(self, state, time, environment, rocketCG):
        airspeed = max(AeroParameters.getLocalFrameAirVel(state, environment).length(), 0.0000001)
        redimConst = self.Lref / (2*airspeed)
        # Calculate moment coefficients from damping coefficients
        localFrameAngularVelocity = Vector(*state.angularVelocity)
        zMomentCoeff = self.zDampingCoeffs * localFrameAngularVelocity * redimConst
        yMomentCoeff = self.yDampingCoeffs * localFrameAngularVelocity * redimConst
        xMomentCoeff = self.xDampingCoeffs * localFrameAngularVelocity * redimConst
        momentCoeffs = [ xMomentCoeff, yMomentCoeff, zMomentCoeff ]

        return AeroFunctions.forceFromCoefficients(state, environment, 0, 0, *momentCoeffs, self.position, self.Aref, self.Lref)

Ancestors

Class variables

var position

Methods

def getAppliedForce(self, state, time, environment, rocketCG)

Expand source code

def getAppliedForce(self, state, time, environment, rocketCG):
    airspeed = max(AeroParameters.getLocalFrameAirVel(state, environment).length(), 0.0000001)
    redimConst = self.Lref / (2*airspeed)
    # Calculate moment coefficients from damping coefficients
    localFrameAngularVelocity = Vector(*state.angularVelocity)
    zMomentCoeff = self.zDampingCoeffs * localFrameAngularVelocity * redimConst
    yMomentCoeff = self.yDampingCoeffs * localFrameAngularVelocity * redimConst
    xMomentCoeff = self.xDampingCoeffs * localFrameAngularVelocity * redimConst
    momentCoeffs = [ xMomentCoeff, yMomentCoeff, zMomentCoeff ]

    return AeroFunctions.forceFromCoefficients(state, environment, 0, 0, *momentCoeffs, self.position, self.Aref, self.Lref)

class AeroForce (componentDictReader, rocket, stage)

A zero-Inertia component with constant aerodynamic coefficients

Expand source code

class AeroForce(RocketComponent):
    ''' A zero-Inertia component with constant aerodynamic coefficients '''
    # Object is just a force, inertia is zero
    inertia = Inertia(Vector(0,0,0), Vector(0,0,0), 0)

    def __init__(self, componentDictReader, rocket, stage):
        self.componentDictReader = componentDictReader
        self.rocket = rocket
        self.stage = stage
        self.name = componentDictReader.getDictName()

        self.position = componentDictReader.getVector("position")
        self.Aref = componentDictReader.getFloat("Aref")
        self.Lref = componentDictReader.getFloat("Lref")

        Cd = componentDictReader.getFloat("Cd")
        Cl = componentDictReader.getFloat("Cl")
        momentCoeffs = componentDictReader.getVector("momentCoeffs")

        self.aeroCoeffs = [ Cd, Cl, *momentCoeffs ]

    def getInertia(self, time, state):
        return self.inertia

    def getAppliedForce(self, state, time, environment, rocketCG):
        return AeroFunctions.forceFromCoefficients(state, environment, *self.aeroCoeffs, self.position, self.Aref, self.Lref)

Ancestors

RocketComponent
abc.ABC

Subclasses

Class variables

var inertia

Methods

def getAppliedForce(self, state, time, environment, rocketCG)

Expand source code

def getAppliedForce(self, state, time, environment, rocketCG):
    return AeroFunctions.forceFromCoefficients(state, environment, *self.aeroCoeffs, self.position, self.Aref, self.Lref)

def getInertia(self, time, state)

Expand source code

def getInertia(self, time, state):
    return self.inertia

class BoatTail (*args)

Defines a conical boattail (aerodynamic properties quite similar to curved boattails, especially in supersonic flight) Always assumes it's at the bottom of a rocket. Modelled like a Transition object, but accounts for base drag.

Two possible sets of inputs:
1. Initialization as a regular, dictionary-defined rocket component:
* args = (componentDictReader, rocket, stage)
* Expected classes: (SubDictReader, Rocket, Stage)
2. Manual initialization:
* args = (startDiameter, endDiameter, length, position, inertia, rocket, stage, name, surfaceRoughness)
* Expected classes: (float, float, float, Vector, Inertia, Rocket, Stage, str, float)

Expand source code

class BoatTail(Transition):
    '''        
        Defines a conical boattail (aerodynamic properties quite similar to curved boattails, especially in supersonic flight)
        Always assumes it's at the bottom of a rocket.
        Modelled like a Transition object, but accounts for base drag.
    '''
    canConnectToComponentBelow = False 
    ''' Overrides attribute inherited from BodyComponent (through Transition), to indicate that this component must exist at the very bottom of a rocket '''

    def getAppliedForce(self, rocketState, time, environment, CG) -> ForceMomentSystem:
        Mach = AeroParameters.getMachNumber(rocketState, environment)
        Aref = self.rocket.Aref
        
        #### Normal Force ####
        AOA = AeroParameters.getTotalAOA(rocketState, environment)
        CN = AeroFunctions.Barrowman_GetCN(AOA, Aref, self.topArea, self.bottomArea)

        #### Pressure Drag ####
        Cd_base = AeroFunctions.getBaseDragCoefficient(Mach)
        Cd_pressure = Cd_base * self.CdAdjustmentFactor
        Cd_pressure *= self.frontalArea / self.rocket.Aref

        noEngine = (self.stage.engineShutOffTime == None)
        if noEngine or time > self.stage.engineShutOffTime:
            # Add base drag if engine is off
            Cd_pressure += Cd_base * self.bottomArea / Aref

        #### Skin Friction Drag ####
        if self.wettedArea == 0:
            skinFrictionDragCoefficient = 0
            rollDampingMoment = Vector(0,0,0)
        else:
            skinFrictionDragCoefficient, rollDampingMoment = AeroFunctions.getCylindricalSkinFrictionDragCoefficientAndRollDampingMoment(rocketState, environment, \
                 self.length, Mach, self.surfaceRoughness, self.wettedArea, Aref, self.rocket.finenessRatio, self.rocket.fullyTurbulentBL)

        #### Total Drag ####
        Cd = Cd_pressure + skinFrictionDragCoefficient

        #### Assemble & return final force object ####
        return AeroFunctions.forceFromCdCN(rocketState, environment, Cd, CN, self.CPLocation, Aref, moment=rollDampingMoment)

Ancestors

Class variables

var canConnectToComponentBelow: Overrides attribute inherited from BodyComponent (through Transition), to indicate that this component must exist at the very bottom of a rocket

Methods

def getAppliedForce(self, rocketState, time, environment, CG) ‑> ForceMomentSystem

Expand source code

def getAppliedForce(self, rocketState, time, environment, CG) -> ForceMomentSystem:
    Mach = AeroParameters.getMachNumber(rocketState, environment)
    Aref = self.rocket.Aref
    
    #### Normal Force ####
    AOA = AeroParameters.getTotalAOA(rocketState, environment)
    CN = AeroFunctions.Barrowman_GetCN(AOA, Aref, self.topArea, self.bottomArea)

    #### Pressure Drag ####
    Cd_base = AeroFunctions.getBaseDragCoefficient(Mach)
    Cd_pressure = Cd_base * self.CdAdjustmentFactor
    Cd_pressure *= self.frontalArea / self.rocket.Aref

    noEngine = (self.stage.engineShutOffTime == None)
    if noEngine or time > self.stage.engineShutOffTime:
        # Add base drag if engine is off
        Cd_pressure += Cd_base * self.bottomArea / Aref

    #### Skin Friction Drag ####
    if self.wettedArea == 0:
        skinFrictionDragCoefficient = 0
        rollDampingMoment = Vector(0,0,0)
    else:
        skinFrictionDragCoefficient, rollDampingMoment = AeroFunctions.getCylindricalSkinFrictionDragCoefficientAndRollDampingMoment(rocketState, environment, \
             self.length, Mach, self.surfaceRoughness, self.wettedArea, Aref, self.rocket.finenessRatio, self.rocket.fullyTurbulentBL)

    #### Total Drag ####
    Cd = Cd_pressure + skinFrictionDragCoefficient

    #### Assemble & return final force object ####
    return AeroFunctions.forceFromCdCN(rocketState, environment, Cd, CN, self.CPLocation, Aref, moment=rollDampingMoment)

Inherited members

Transition:
- getBottomInterfaceLocation
- getMaxDiameter
- getRadius
- getTopInterfaceLocation

class BodyComponent

Class that defines interface for axisymmetric body components. Contains logic for detecting adjacent body components & defining interfaces with them Examples: NoseCone, Stage

Expand source code

class BodyComponent(ABC):
    ''' 
        Class that defines interface for axisymmetric body components.
        Contains logic for detecting adjacent body components & defining interfaces with them 
        Examples: `MAPLEAF.Rocket.NoseCone`, `MAPLEAF.Rocket.Stage`
    '''
    # Override these attributes in child classes to change whether they can connect to components above/below them
    canConnectToComponentAbove = True
    canConnectToComponentBelow = True

    def getTopInterfaceLocation(self) -> Union[None, Vector]:
        ''' For planar cylindrical interfaces, returns the location of the center of the cylindrical interface '''
        if self.canConnectToComponentAbove:
            return self.position
        else:
            return None

    def getBottomInterfaceLocation(self) -> Union[None, Vector]:
        ''' For planar cylindrical interfaces, returns the location of the center of the cylindrical interface '''
        if self.canConnectToComponentBelow:
            baseZCoord = self.position.Z-self.length
            return Vector(self.position.X, self.position.Y, baseZCoord)
        else:
            return None

    def _getCenterOfPressure(self, *args) -> Vector:
        return self.CPLocation

    @abstractmethod
    def getMaxDiameter(self):
        ''' These functions used for determining rocket's current max diameter '''
        return

    @abstractmethod
    def getRadius(self, distanceFromTop: float) -> float:
        ''' Should return body component radius as a function of distance from the top of the component '''
        return

Ancestors

abc.ABC

Subclasses

Class variables

var canConnectToComponentAbove
var canConnectToComponentBelow

Methods

def getBottomInterfaceLocation(self) ‑> Optional[None]

For planar cylindrical interfaces, returns the location of the center of the cylindrical interface

Expand source code

def getBottomInterfaceLocation(self) -> Union[None, Vector]:
    ''' For planar cylindrical interfaces, returns the location of the center of the cylindrical interface '''
    if self.canConnectToComponentBelow:
        baseZCoord = self.position.Z-self.length
        return Vector(self.position.X, self.position.Y, baseZCoord)
    else:
        return None

def getMaxDiameter(self)

These functions used for determining rocket's current max diameter

Expand source code

@abstractmethod
def getMaxDiameter(self):
    ''' These functions used for determining rocket's current max diameter '''
    return

def getRadius(self, distanceFromTop: float) ‑> float

Should return body component radius as a function of distance from the top of the component

Expand source code

@abstractmethod
def getRadius(self, distanceFromTop: float) -> float:
    ''' Should return body component radius as a function of distance from the top of the component '''
    return

def getTopInterfaceLocation(self) ‑> Optional[None]

For planar cylindrical interfaces, returns the location of the center of the cylindrical interface

Expand source code

def getTopInterfaceLocation(self) -> Union[None, Vector]:
    ''' For planar cylindrical interfaces, returns the location of the center of the cylindrical interface '''
    if self.canConnectToComponentAbove:
        return self.position
    else:
        return None

class BodyTube (componentDictReader, rocket, stage)

Represent a cylindrical body tube

Expand source code

class BodyTube(FixedMass, BodyComponent):
    ''' Represent a cylindrical body tube '''
    
    #### Init Functions ####
    def __init__(self, componentDictReader, rocket, stage):
        FixedMass.__init__(self, componentDictReader, rocket, stage)

        self.length = componentDictReader.getFloat("length")
        self.outerDiameter = componentDictReader.getFloat("outerDiameter")
        self.surfaceRoughness = componentDictReader.tryGetFloat("surfaceRoughness", defaultValue=self.rocket.surfaceRoughness)

        # Tell the rocket/stage what its diameter is
        self.rocket.bodyTubeDiameter = self.outerDiameter
        self.stage.bodyTubeDiameter = self.outerDiameter

        self._precomputeGeometry()

    def _precomputeGeometry(self):
        self.volume = math.pi * self.outerDiameter**2 / 4 * self.length
        self.planformArea = self.outerDiameter * self.length
        self.wettedArea = math.pi * self.outerDiameter * self.length
        
        self.CPLocation = self.position - Vector(0,0,self.length/2)

    #### Operational Functions ####
    def getAppliedForce(self, rocketState, time, environment, CG):
        Aref = self.rocket.Aref

        # Normal Force ----------------------------------------------------------------------------------------
        AOA = AeroParameters.getTotalAOA(rocketState, environment)
        # Niskanen Eqn 3.26 - originally from Galejs
        normalForceCoefficient = 1.1 * (math.sin(AOA))**2
        normalForceCoefficient *= self.planformArea / Aref

        # Drag -----------------------------------------------------------------------------------------------
        # Skin Friction
        Mach = AeroParameters.getMachNumber(rocketState, environment)
        skinFrictionDragCoefficient, rollDampingMoment = AeroFunctions.getCylindricalSkinFrictionDragCoefficientAndRollDampingMoment(rocketState, environment, self.length, Mach, self.surfaceRoughness, \
            self.wettedArea, Aref, self.rocket.finenessRatio, self.rocket.fullyTurbulentBL)
        
        #There is no pressure drag associated with bodytubes - skin friction drag is total drag

        # Damping moments --------------------------------------------------------------------------------------
        dampingMoments = self._computeLongitudinalDampingMoments(rocketState, environment, CG)
        
        # Roll damping due to skin friction
        dampingMoments += rollDampingMoment

        # Compute & dimensionalize total-------------------------------------------------------------------------
        return AeroFunctions.forceFromCdCN(rocketState, environment, skinFrictionDragCoefficient, normalForceCoefficient, self.CPLocation, Aref, moment=dampingMoments)

    def _computeLongitudinalDampingMoments(self, rocketState, environment, CoR, nSegments=25):
        ''' CoR = Center of Rotation (usually the rocket's CG) '''
        # TODO: Check if there's an nice analytical solution to this integral

        # Eqn (3.57) from Niskanen, converted to vector form and integrated numerically
        # Numerical integration: Divide body tube into segments
            # Works for body tubes in any position, rotating on any axis
        dL = self.length / nSegments
        dA = dL * self.outerDiameter
        totalDampingMoment = Vector(0,0,0)

        for i in range(nSegments):
            # Center of current segment of body tube
            segmentCentroid = self.position + Vector(0,0, -dL*i - dL/2)
            
            # Vector from center of rotation to segment centroid
            psi = segmentCentroid - CoR 

            # Velocity due to rotation is angular velocity cross distance (from center of rotation) vector
            segmentVelocity = rocketState.angularVelocity.crossProduct(psi)
            
            # Re-dimensionalizing coefficient required squared velocity - while preserving sign
            sqVel = Vector( segmentVelocity.X*abs(segmentVelocity.X), segmentVelocity.Y*abs(segmentVelocity.Y), \
                 segmentVelocity.Z*abs(segmentVelocity.Z) )

            # Moment is distance vector cross force vector
            dM = -psi.crossProduct(sqVel)
            totalDampingMoment += dM

        # All terms of the integral summation must be multiplied by this constant:
            # 1.1 is drag coefficient of a cylinder in crossflow
            # 1/2 and density are from re-dimensionalizing the coefficient
            # Outerdiameter is part of area calculation
        integrationConstantMultiple = (1.1/2)*environment.Density*dA
        totalDampingMoment *= integrationConstantMultiple

        return totalDampingMoment

    def getMaxDiameter(self):
        return self.outerDiameter

    def getRadius(self, _):
        return self.outerDiameter/2

    def plotShape(self):
        import matplotlib.pyplot as plt
        
        # Assume body tube is on the Z-axis
        tipZ = self.position.Z
        tailZ = tipZ - self.length
        radius = self.outerDiameter/2

        Xvals = []
        Yvals = []
        Xvals.append(tipZ)
        Yvals.append(radius) # top right

        Xvals.append(tipZ)
        Yvals.append(-radius) # bottom right

        Xvals.append(tailZ)
        Yvals.append(-radius) # bottom left 

        Xvals.append(tailZ)
        Yvals.append(radius) # top left

        Xvals.append(tipZ)
        Yvals.append(radius) # close the square
        plt.plot(Xvals, Yvals, color = 'k')

Ancestors

Methods

def getAppliedForce(self, rocketState, time, environment, CG)

Expand source code

def getAppliedForce(self, rocketState, time, environment, CG):
    Aref = self.rocket.Aref

    # Normal Force ----------------------------------------------------------------------------------------
    AOA = AeroParameters.getTotalAOA(rocketState, environment)
    # Niskanen Eqn 3.26 - originally from Galejs
    normalForceCoefficient = 1.1 * (math.sin(AOA))**2
    normalForceCoefficient *= self.planformArea / Aref

    # Drag -----------------------------------------------------------------------------------------------
    # Skin Friction
    Mach = AeroParameters.getMachNumber(rocketState, environment)
    skinFrictionDragCoefficient, rollDampingMoment = AeroFunctions.getCylindricalSkinFrictionDragCoefficientAndRollDampingMoment(rocketState, environment, self.length, Mach, self.surfaceRoughness, \
        self.wettedArea, Aref, self.rocket.finenessRatio, self.rocket.fullyTurbulentBL)
    
    #There is no pressure drag associated with bodytubes - skin friction drag is total drag

    # Damping moments --------------------------------------------------------------------------------------
    dampingMoments = self._computeLongitudinalDampingMoments(rocketState, environment, CG)
    
    # Roll damping due to skin friction
    dampingMoments += rollDampingMoment

    # Compute & dimensionalize total-------------------------------------------------------------------------
    return AeroFunctions.forceFromCdCN(rocketState, environment, skinFrictionDragCoefficient, normalForceCoefficient, self.CPLocation, Aref, moment=dampingMoments)

def plotShape(self)

Expand source code

def plotShape(self):
    import matplotlib.pyplot as plt
    
    # Assume body tube is on the Z-axis
    tipZ = self.position.Z
    tailZ = tipZ - self.length
    radius = self.outerDiameter/2

    Xvals = []
    Yvals = []
    Xvals.append(tipZ)
    Yvals.append(radius) # top right

    Xvals.append(tipZ)
    Yvals.append(-radius) # bottom right

    Xvals.append(tailZ)
    Yvals.append(-radius) # bottom left 

    Xvals.append(tailZ)
    Yvals.append(radius) # top left

    Xvals.append(tipZ)
    Yvals.append(radius) # close the square
    plt.plot(Xvals, Yvals, color = 'k')

Inherited members

BodyComponent:
- getBottomInterfaceLocation
- getMaxDiameter
- getRadius
- getTopInterfaceLocation

class EventTypes (value, names=None, *, module=None, qualname=None, type=None, start=1)

An enumeration.

Expand source code

class EventTypes(Enum):
    Apogee = "apogee"
    AscendingThroughAltitude = "ascendingThroughAltitude"
    DescendingThroughAltitude = "descendingThroughAltitude"
    MotorBurnout = "motorBurnout"
    TimeReached = "timeReached"

Ancestors

enum.Enum

Class variables

var Apogee
var AscendingThroughAltitude
var DescendingThroughAltitude
var MotorBurnout
var TimeReached

class Fin (componentDictReader, parentFinSet, spanwiseDirection, rocket, stage)

Represents a single fin in a FinSet

Expand source code

class Fin(FixedMass):
    ''' Represents a single fin in a FinSet '''

    def __init__(self, componentDictReader, parentFinSet, spanwiseDirection, rocket, stage):

        FixedMass.__init__(self, componentDictReader, rocket, stage)

        self.finSet = parentFinSet # type: FinSet
        ''' Reference to the parent fin set '''

        self.spanSliceAreas = self.finSet.spanSliceAreas
        self.spanSliceRadii = self.finSet.spanSliceRadii
        self.sliceLEPositions = self.finSet.sliceLEPositions
        self.sliceLengths = self.finSet.sliceLengths

        self.finAngle = self.finSet.finAngle
        self.stallAngle = self.finSet.stallAngle

        self.span = self.finSet.span
        self.planformArea = self.finSet.planformArea
        self.midChordSweep = self.finSet.midChordSweep
        
        # Unit vector parallel to the fin's shaft and/or fin's spanwise direction, normal to the rocket's body tube
        self.spanwiseDirection = spanwiseDirection

        # Vector normal to the fin surface, in the X/Y Plane, when self.finAngle == 0
        self.undeflectedFinNormal = Vector(self.spanwiseDirection[1], -self.spanwiseDirection[0], 0)

        self.tipPosition = self.finSet.dZ_dSpan_LE_neg * self.span
        
        # Spanwise component of Center of Pressure (CP) location
        self.CPSpanWisePosition = self.spanwiseDirection*(self.finSet.MACYPos + self.finSet.bodyRadius)

    def _barrowmanAeroFunc(self, rocketState, time, environment, precomputedData, CG = Vector(0,0,0), finDeflectionAngle = None):
        '''
            Precomputed Data is a named tuple (PreComputedFinAeroData) which contains data/results from the parts of the fin aerodynamics calculation
                that are common to all fins in a FinSet (drag calculations mostly). These calculations are performed at the FinSet level.
            Only the parts of the Fin Aero Computation that change from fin to fin (normal force mostly, since different fins can have different angles of attack) are computed here
        '''
        Mach = AeroParameters.getMachNumber(rocketState, environment)
        dynamicPressure = AeroParameters.getDynamicPressure(rocketState, environment)

        if finDeflectionAngle == None:
            finDeflectionAngle = self.finAngle # Adjusted by parent finset during each timestep, when the FinSet is controlled

        # Unpack precomputedData
        airVelRelativeToFin, CPXPos, totalDragCoefficient = precomputedData

        #### Compute normal force-----------------------------------------------------------------------------------------------------------------
           
        # Find fin normal vector after fin deflection       
        finDeflectionRotation = Quaternion(axisOfRotation=self.spanwiseDirection, angle=radians(finDeflectionAngle))
        finNormal = finDeflectionRotation.rotate(self.undeflectedFinNormal)
        # Get the tangential velocity component vector, per unit radial distance from the rocket centerline
        rollAngVel = AngularVelocity(0, 0, rocketState.angularVelocity.Z)
        unitSpanTangentialAirVelocity = rollAngVel.crossProduct(self.spanwiseDirection)*(-1)
        
        def subsonicNormalForce(Mach): # Subsonic linear method
            tempBeta = AeroParameters.getBeta(Mach)
            CnAlpha = getFinCnAlpha_Subsonic_Barrowman(self.span, self.planformArea, tempBeta, self.midChordSweep)
            return getSubsonicFinNormalForce(airVelRelativeToFin, unitSpanTangentialAirVelocity, finNormal, self.spanwiseDirection, self.CPSpanWisePosition.length(), CnAlpha, self)

        def supersonicNormalForce(Mach): # Supersonic Busemann method
            gamma = AeroFunctions.getGamma()
            tempBeta = AeroParameters.getBeta(Mach)
            K1, K2, K3, Kstar = getBusemannCoefficients(Mach, tempBeta, gamma)

            # Mach Cone coords
            machAngle = asin(1/Mach)
            machCone_negZPerRadius = 1 / tan(machAngle)
            machConeEdgeZPos = []
            outerRadius = self.spanSliceRadii[-1]
            for i in range(len(self.spanSliceRadii)):
                machConeAtCurrentRadius = (outerRadius - self.spanSliceRadii[i])*machCone_negZPerRadius + self.tipPosition
                machConeEdgeZPos.append(machConeAtCurrentRadius)

            return getSupersonicFinNormalForce(airVelRelativeToFin, unitSpanTangentialAirVelocity, finNormal, machConeEdgeZPos, self.spanwiseDirection, self.CPSpanWisePosition.length(), K1, K2, K3, Kstar, self)

        if Mach <= 0.8:
            normalForceMagnitude, finMoment = subsonicNormalForce(Mach)

        elif Mach < 1.4:
            # Interpolate in transonic region
            # TODO: Do this with less function evaluations? Perhaps precompute AOA and Mach combinations and simply interpolate? Lazy precompute? Cython?
            x1, x2 = 0.8, 1.4 # Start, end pts of interpolated region
            dx = 0.001

            # Find normal force and derivative at start of interpolation interval
            f_x1, m_x1 = subsonicNormalForce(x1)
            f_x12, m_x12 = subsonicNormalForce(x1+dx)

            # Find normal force and derivative at end of interpolation interval
            f_x2, m_x2 = supersonicNormalForce(x2)
            f_x22, m_x22 = supersonicNormalForce(x2+dx)

            normalForceMagnitude = Interpolation.cubicInterp(Mach, x1, x2, f_x1, f_x2, f_x12, f_x22, dx)
            finMoment = Interpolation.cubicInterp(Mach, x1, x2, m_x1, m_x2, m_x12, m_x22, dx)
        else:
            normalForceMagnitude, finMoment = supersonicNormalForce(Mach)
             
        # Complete redimensionalization of normal force coefficients by multiplying by dynamic pressure
        # Direct the normal force along the fin's normal direction
        normalForce = normalForceMagnitude * dynamicPressure * finNormal
        finMoment *= dynamicPressure
                
        #### Get axial force-----------------------------------------------------------------------------------------------------------------------
        
        avgAOA = getFinSliceAngleOfAttack(self.spanSliceRadii[round(len(self.spanSliceAreas)/2)], airVelRelativeToFin, unitSpanTangentialAirVelocity, finNormal, self.spanwiseDirection, self.stallAngle) # Approximate average fin AOA
        totalAxialForceCoefficient = AeroFunctions.getDragToAxialForceFactor(avgAOA) * totalDragCoefficient
        axialForceMagnitude = totalAxialForceCoefficient * self.rocket.Aref * dynamicPressure
        axialForceDirection = self.spanwiseDirection.crossProduct(finNormal)
        axialForce = axialForceDirection*axialForceMagnitude
        
        #### Get CP Location ----------------------------------------------------------------------------------------------------------------------

        CPChordWisePosition = self.position - Vector(0, 0, CPXPos) # Ignoring the change in CP position due to fin deflection for now
        globalCP = self.CPSpanWisePosition + CPChordWisePosition

        #### Assemble total force moment system objects--------------------------------------------------------------------------------------------

        totalForce = normalForce + axialForce
        return ForceMomentSystem(totalForce, globalCP, moment=Vector(0, 0, finMoment)), globalCP

    def getAppliedForce(self, rocketState, time, environment, CG, precomputedData):
        [ aeroForce, CP ] = self._barrowmanAeroFunc(rocketState, time, environment, precomputedData, CG)
        return aeroForce

Ancestors

Instance variables

var finSet: Reference to the parent fin set

Methods

def getAppliedForce(self, rocketState, time, environment, CG, precomputedData)

Expand source code

def getAppliedForce(self, rocketState, time, environment, CG, precomputedData):
    [ aeroForce, CP ] = self._barrowmanAeroFunc(rocketState, time, environment, precomputedData, CG)
    return aeroForce

class FinSet (componentDictReader, rocket, stage)

Class represents a set of n identical fins, all at the same longitudinal location, arranged axisymmetrically. Fin orientations can be controlled by a ControlSystem

Expand source code

class FinSet(FixedMass, ActuatedSystem):
    ''' Class represents a set of n identical fins, all at the same longitudinal location, arranged axisymmetrically. Fin orientations can be controlled by a `MAPLEAF.GNC.ControlSystems.ControlSystem` '''

    #### Initialization Functions ####
    def __init__(self, componentDictReader, rocket, stage):
        FixedMass.__init__(self, componentDictReader, rocket, stage)

        self.controlSystem = None   # Set by self.initializeActuators, if it is called by a control system being initialized ex: RocketControlSystem.__init__()
        self.actuatorList = None    # Set by self.initializeActuators, if it is called by a control system being initialized ex: RocketControlSystem.__init__()

        self.numFins = componentDictReader.getInt("numFins")
        self.firstFinAngle = componentDictReader.getFloat("firstFinAngle")

        # Set fin properties
        self.finAngle = componentDictReader.getFloat("finCantAngle")
        self.sweepAngle = math.radians(componentDictReader.getFloat("sweepAngle"))
        self.rootChord = componentDictReader.getFloat("rootChord")
        self.tipChord = componentDictReader.getFloat("tipChord")
        self.span = componentDictReader.getFloat("span")
        self.surfaceRoughness = componentDictReader.tryGetFloat("surfaceRoughness", defaultValue=self.rocket.surfaceRoughness)
        self.thickness = componentDictReader.getFloat("thickness")
        self.numFinSpanSteps = componentDictReader.getInt("numFinSpanSlicesForIntegration")
        
        # Leading Edge properties
        self.leadingEdgeShape = componentDictReader.getString("LeadingEdge.shape")
        if self.leadingEdgeShape == "Round":
            self.leadingEdgeRadius = componentDictReader.tryGetFloat("LeadingEdge.radius", defaultValue=self.thickness/100)
        elif self.leadingEdgeShape == "Blunt":
            self.leadingEdgeThickness = componentDictReader.tryGetFloat("LeadingEdge.thickness", defaultValue=self.thickness)
        else:
            raise ValueError("ERROR: Leading edge shape: {} not implemented. Use 'Round' or 'Blunt'".format(self.leadingEdgeShape))

        # Trailing Edge properties
        self.trailingEdgeShape = componentDictReader.getString("TrailingEdge.shape")
        if self.trailingEdgeShape == "Round":
            self.trailingEdgeRadius = componentDictReader.getFloat("TrailingEdge.radius")
        elif self.trailingEdgeShape == "Blunt":
            self.trailingEdgeThickness = componentDictReader.getFloat("TrailingEdge.thickness")
        elif self.trailingEdgeShape != "Tapered":
            raise ValueError("ERROR: Trailing edge shape: {} not implemented. Use 'Round', 'Blunt', or 'Tapered'".format(self.trailingEdgeShape))
        

    def initializeActuators(self, controlSystem):
        self.controlSystem = controlSystem

        # Initialize an actuator model for each fin
        ActuatedSystem.__init__(self, self.numFins)

    def precomputeProperties(self):
        self._calculateSweepAngles()

        ### Compute Other Simple Properties ####
        self.planformArea = self.span * (self.tipChord + self.rootChord) / 2
        self.wettedArea = 2*self.planformArea
        self.aspectRatio = (2*self.span)**2/self.planformArea # Aspect ratio of the wing that would be created by reflecting the fin about its root chord (Fleeman/Niskanen)
        self.stallAngle = 15 * (1 + 1.1/self.aspectRatio) # Compute stall angle from Hoerner, 'Fluid Dynamic Lift' for low aspect ratio wings

        # Get the body radius where the fin is mounted
        stageTopZ = self.stage.position.Z
        distanceFromTopOfStageToTopOfThisFinSet = stageTopZ - self.position.Z
        self.bodyRadius = self.stage.getRadius(distanceFromTopOfStageToTopOfThisFinSet)

        self._calculateInterferenceFactors()
        self._calculateMACProperties()
        self._splitFinsIntoSlices()

        self._precomputeSubsonicFinThicknessDrag()
        self._precomputeCPInterpolationPolynomial()

        # Initialize all the child fins in the self.finList list
        self._initChildFins(self.componentDictReader, self.rocket, self.stage)
    
    def _calculateSweepAngles(self):
        ''' Compute Trailing Edge (TE) and Mid Chord sweep angles '''
        LEtipChordBehindRootChord = self.span * math.tan(self.sweepAngle) # Z-distance the front of the tip chord is behind the front of the root chord
        TETipChordBehindRootChord = self.tipChord + LEtipChordBehindRootChord - self.rootChord # Z-distance the back of the tip chord is behind the back of the root chord

        if(TETipChordBehindRootChord == 0):
            self.trailingEdgeSweep = 0.0
        elif(TETipChordBehindRootChord < 0):
            # Negative trailing edge sweep is a forward-swept rear
            self.trailingEdgeSweep = -1*math.atan2(abs(TETipChordBehindRootChord),self.span)
        else:
            self.trailingEdgeSweep = math.atan2(abs(TETipChordBehindRootChord), self.span)

        self.midChordSweep = (self.sweepAngle + self.trailingEdgeSweep) / 2

        self.dZ_dSpan_LE_neg = math.tan(self.sweepAngle) # -dZ/dSpan (leading edge slope)
        self.dZ_dSpan_TE_neg = math.tan(self.trailingEdgeSweep) # dZ/dSpan (trailing edge slope)

    def _calculateInterferenceFactors(self):
        # Interference w/ body tube
        # Accounts for additional normal forces generated by the fins, due to the presence of the rocket body
        # Niskanen Eqn 3.56, originally from Barrowman
        self.bodyOnFinInterferenceFactor = 1 + self.bodyRadius / (self.span + self.bodyRadius)

        # Interference b/w the fins
        # Niskanen Eqn 3.54, from MIL-HDBK-762
        # Accounts for reduced normal force generated by large numbers of fins in a group
        # Not expected to work well for more than 8 fins
        if self.numFins <= 4:
            self.finNumberInterferenceFactor = 1.0
        elif self.numFins <= 8:
            numberToFactorMap = {
                5: 0.948,
                6: 0.913,
                7: 0.854,
                8: 0.810,
            }
            self.finNumberInterferenceFactor = numberToFactorMap[self.finNumber]
        else:
            self.finNumberInterferenceFactor = 0.75

    def _calculateMACProperties(self):
        ''' Calculates MAC Length, MACY, and MACX - Niskanen (Eqn. 3.30-3.32) '''
        def LEdgeAtY(y): # Leading Edge at y - where y is the spanwise coordinate
            return self.dZ_dSpan_LE_neg*y
                    
        numSpanSteps = round(self.span / 0.0001) # One step per 0.1 mm
        finSpanStep = self.span / numSpanSteps

        MACIntegralSum = 0
        MACYPosIntegralSum = 0
        XMACLeadingEdgeIntegralSum = 0
        for i in range(numSpanSteps):
            y = i * finSpanStep + finSpanStep/2
            MACIntegralSum += ((self.getChord(y)**2)*finSpanStep) # Eqn 3.30, integrated numerically
            MACYPosIntegralSum += y * (self.getChord(y) * finSpanStep) # Eqn 3.31, integrated numerically
            XMACLeadingEdgeIntegralSum += LEdgeAtY(y) * self.getChord(y)*finSpanStep # Eqn 3.32, integrated numerically

        # Save results
        self.MACLength = MACIntegralSum/self.planformArea
        self.MACYPos = MACYPosIntegralSum/self.planformArea
        self.XMACLeadingEdge = XMACLeadingEdgeIntegralSum/self.planformArea

    def _splitFinsIntoSlices(self):
        ''' Precomputes fin area slices for normal force integration '''
        stepSize = self.span / self.numFinSpanSteps
        bodyRadius = self.bodyRadius
        
        self.spanSliceAreas = [] # Area of each slice, m^2
        self.spanSliceRadii = [] # Radius of the center of each slice, from the rocket centerline, m
        self.sliceLEPositions = [] # Z - Position of the LE at the center of the slice, from the tip of root chord
        self.sliceLengths = []

        for i in range(self.numFinSpanSteps):
            y = i*stepSize + stepSize*0.5

            self.spanSliceRadii.append(y + bodyRadius)
            self.spanSliceAreas.append(stepSize*self.getChord(y))
            self.sliceLEPositions.append(y * self.dZ_dSpan_LE_neg)
            self.sliceLengths.append(self.getChord(y))

    def _precomputeSubsonicFinThicknessDrag(self):
        ''' Precompute the subsonic thickness drag coefficient (Barrowman Eqn 4-36) '''
        sigma = (self.aspectRatio/2) * (self.thickness / self.rootChord)**(1/3) # Barrowman Eqn 4-35 b
        CD_TT_star = 1.15 * (self.thickness/self.rootChord)**(5/3) * (1.61 + sigma - sqrt((sigma - 1.43)**2 + 0.578))
        
        # Get skin friction coefficient at Mach 1
        mach1State = copy.deepcopy(self.rocket.rigidBody.state)
        mach1State.velocity = Vector(0,0,340.3) # Get speed up around Mach 1
        environment = self.rocket.environment.getAirProperties(self.rocket.rigidBody.state.position, 0)
        CD_f_star = AeroFunctions.getSkinFrictionCoefficient(mach1State, environment, self.MACLength, 1.0, self.surfaceRoughness, self.rocket.fullyTurbulentBL)
        
        # Final result
        self.subsonicFinThicknessK = cos(self.midChordSweep)**2 + (((CD_TT_star/CD_f_star - 4*cos(self.midChordSweep)*(self.thickness/self.rootChord)) / (120 * cos(self.midChordSweep)**2 * (self.thickness/self.rootChord)**4))**2)**(1/3)

    def _precomputeCPInterpolationPolynomial(self):
        ''' Precompute Interpolation polynomial coefficients for getCPXLocation function '''
        AR = self.aspectRatio
        f_2 = (AR * (3)**0.5 - 0.67) / (2 * (3)**0.5 * AR - 1)
        f_prime_2_t1 = (2 * (3)**-0.5 * AR) * (2 * AR * (3)**0.5 - 1) 
        f_prime_2_t2 = ((3)**0.5 * AR - 0.67) * (AR * 4 * (3)**-0.5)
        f_prime_2_denom = (2 * AR * (3)**0.5 -1)**2
        f_prime_2 = (f_prime_2_t1 - f_prime_2_t2) / f_prime_2_denom
        
        a_mat = [
                [0.5**5, 0.5**4, 0.5**3, 0.5**2, 0.5, 1],
                [5*(0.5)**4, 4*(0.5)**3, 3*(0.5)**2, 2*(0.5), 1, 0],
                [2**5, 2**4, 2**3, 2**2, 2, 1],
                [5*(2)**4, 4*(2)**3, 3*(2)**2, 2*2, 1, 0],
                [20*2**3, 12*2**2, 6*2, 2, 0, 0],
                [60*2**2, 24*2, 6, 0, 0, 0],
                ]
        
        b_vec = [0.25, 0, f_2, f_prime_2, 0, 0]
        
        # Save results
        self.x = np.linalg.solve(a_mat, b_vec)   

    def _initChildFins(self, componentDictReader, rocket, stage):
        self.finList = []
        
        finSeparationAngle = 360/self.numFins
        firstFinAngle = self.firstFinAngle

        # Create child fins
        for i in range(self.numFins):
            # Figure out what direction the fin points in
            spanwiseUnitX = math.cos(math.radians(firstFinAngle + i*finSeparationAngle))
            spanwiseUnitY = math.sin(math.radians(firstFinAngle + i*finSeparationAngle))
            spanWiseDirection = Vector(spanwiseUnitX, spanwiseUnitY, 0)

            #Initialize each fin separately and keep in a list
            self.finList.append(Fin(componentDictReader, self, spanWiseDirection, rocket, stage))

    # TODO: Replace this with control system log
    # TODO: Remove all getLogHeader functions
    def getLogHeader(self):
        header = " {}FX(N) {}FY(N) {}FZ(N) {}MX(Nm) {}MY(Nm) {}MZ(Nm)".format(*[self.name]*6)
        
        if self.controlSystem != None:
            for finNumber in range(1, self.numFins+1):
                header += " {}Fin{}Angle(deg)".format(self.name, finNumber)

        return header

    ### Functions used during simulation ###
    def getAppliedForce(self, rocketState, time, environment, CG):
        #### If control system exists, use actuator deflections 1:1 to set fin angles ####
        if self.controlSystem != None:
            # Update fin angles
            for i in range(self.numFins):
                self.finList[i].finAngle = self.actuatorList[i].getDeflection(time)

        #### Pre-calculate common properties for all child Fins ####
        precomputedData = self._getPreComputedFinAeroData(rocketState, environment, CG)

        #### Add up forces from all child Fins ####
        aeroForce = ForceMomentSystem(Vector(0,0,0), self.position)
        for fin in self.finList:
            aeroForce += fin.getAppliedForce(rocketState, time, environment, CG, precomputedData)

        # TODO: Correct for sub/transonic rolling moment fin-fin interference from a high number of fins
        
        # Apply fin-body interference factor to total forces in the normal (X/Y) directions and moments
        totalInterferenceFactor = self.bodyOnFinInterferenceFactor * self.finNumberInterferenceFactor
        aeroForce.force.X *= totalInterferenceFactor
        aeroForce.force.Y *= totalInterferenceFactor
        aeroForce.moment.X *= totalInterferenceFactor
        aeroForce.moment.Y *= totalInterferenceFactor

        return aeroForce

    def _getPreComputedFinAeroData(self, rocketState, environment, CG):
        #General Info ---------------------------------------------------------------------------------------------------------------------
        Aref = self.rocket.Aref
        Mach = AeroParameters.getMachNumber(rocketState, environment)
        dynamicPressure = AeroParameters.getDynamicPressure(rocketState, environment)

        # Skin Friction Drag -------------------------------------------------------------------------------------------------------------
        
        skinFrictionCoefficient = AeroFunctions.getSkinFrictionCoefficient(rocketState, environment, self.MACLength, Mach, self.surfaceRoughness, self.rocket.fullyTurbulentBL)    
        # Adjust to the rocket reference area (skinFrictionCoefficient is based on wetted area)
        skinFrictionDragCoefficient = skinFrictionCoefficient*(self.wettedArea / Aref)
        # Correct for additional surface area due to fin thickness - Niskanen Eqn 3.85
        skinFrictionDragCoefficient *= (1 + 2*self.thickness/self.MACLength)
        
        # Pressure Drag ------------------------------------------------------------------------------------------------------------------

        # Leading edge drag coefficient
        if self.leadingEdgeShape == "Round":
            leadingEdgeCd = AeroFunctions.getCylinderCrossFlowCd_ZeroBaseDrag(Mach)
            LEthickness = self.leadingEdgeRadius*2
        elif self.leadingEdgeShape == "Blunt":
            leadingEdgeCd = AeroFunctions.getBluntBodyCd_ZeroBaseDrag(Mach)
            LEthickness = self.leadingEdgeThickness

        # Adjust for leading edge angle and convert reference area to rocket reference area - Barrowman Eqn 4-22
        leadingEdgeCdAdjustmentFactor = LEthickness * self.span * cos(self.sweepAngle)**2 / Aref
        leadingEdgeCd *= leadingEdgeCdAdjustmentFactor
            
        # Trailing edge drag coefficient - simpler method from Niskanen section 3.4.4. Corrected to use only the trailing edge area, not the full fin frontal area
            # more intricate method available in Barrowman
        baseDragCd = AeroFunctions.getBaseDragCoefficient(Mach)
        if self.trailingEdgeShape == "Tapered":
            TEthickness = 0 # Zero base drag
        elif self.trailingEdgeShape == "Round":
            TEthickness = self.trailingEdgeRadius # 0.5 base drag
        elif self.trailingEdgeShape == "Blunt":
            TEthickness = self.trailingEdgeThickness # Full base drag
        
        # Convert to standard rocket reference Area
        trailingEdgeCd = baseDragCd * self.span*TEthickness/Aref

        # Thickness / Wave Drag
        #TODO: This section doesn't seem to be working quite right
        if Mach <= 1:
            # Thickness drag, subsonic
            thicknessDrag = 4*skinFrictionCoefficient*((self.thickness/self.rootChord)*cos(self.midChordSweep) + \
                (30 * (self.thickness/self.rootChord)**4 * cos(self.midChordSweep)**2) /  \
                    (self.subsonicFinThicknessK - Mach**2 * cos(self.midChordSweep)**2)**(3/2))
        else:
            # Supersonic wave drag
            # Using simplistic method from Hoerner - assumes diamond profile (pg 17-12, Eqn 29)
            # TODO: Implement method from Barrowman's FIN program
            thicknessDrag = 2.3*self.aspectRatio * (self.thickness / self.MACLength)**2

        thicknessDrag *= self.planformArea/Aref

        pressureDragCoefficient = leadingEdgeCd + trailingEdgeCd + thicknessDrag
                        
        # Total Drag --------------------------------------------------------------------------------------------------------------------
        
        totalDragCoefficient = pressureDragCoefficient + skinFrictionDragCoefficient

        localFrameRocketVelocity = AeroParameters.getLocalFrameAirVel(rocketState, environment)
        axialPositionRelCG = self.position - CG
        finVelocityDueToRocketPitchYaw = rocketState.angularVelocity.crossProduct(axialPositionRelCG) 
        airVelRelativeToFin = localFrameRocketVelocity - finVelocityDueToRocketPitchYaw # The negative puts it in the wind frame
        CPXPos = self._getCPXPos(Mach)

        #### Transfer info to fins ####
        return PreComputedFinAeroData(airVelRelativeToFin, CPXPos, totalDragCoefficient)

    def _getCPXPos(self, Mach):
        ''' 
            Finds the chordwise position of the center of pressure, relative to the tip of the root chord (self.position) 
            Methods from Tactical missile design by Fleeman through Niskanen
        '''
        if Mach < 0.5:
            XF = self.XMACLeadingEdge + 0.25*self.MACLength #As per open rocket documentation

        elif Mach >= 2.0:
            beta = AeroParameters.getBeta(Mach)
            CPCalc = lambda AR, B : (AR*B - 0.67) / (2*AR*B - 1) #Eq 3.35 of open rocket documentation (greater than M = 2)
            XF = self.MACLength*CPCalc(self.aspectRatio, beta) + self.XMACLeadingEdge #As per open rocket documentation

        else:
            #Between the two extremes is a polynomial curve fit
            #Fifth order polynomial fit as per Niskanen Eqn 3.36 (Originally from Fleeman)
                # Old code: XF = np.polyval(self.x, Mach)*self.MACLength + self.XMACLeadingEdge
                # Evaluate polynomial manually for maximum speed (8.5x improvement for 6 coeffs)
                    # Speed test in test/speedTests/PolyvalSpeed.py
            polyval = 0
            nTerms = len(self.x)
            for i in range(nTerms):
                polyval += self.x[i]*Mach**(nTerms-1-i)

            XF = polyval*self.MACLength + self.XMACLeadingEdge # Polyval defines the fraction of the MAC that the Cp is behind the leading edge

        return XF

    def getChord(self, y):
        ''' Returns chord at y, where y==0 is at the root chord, and y==span is at the tip chord '''
        return (self.tipChord - self.rootChord) * (y/self.span) + self.rootChord

    def plotShape(self):
        import matplotlib.pyplot as plt

        BodyTubeDiameter = self.bodyRadius*2
        SubComponentPos = self.position.Z
        span = self.span
        rootChord = self.rootChord
        tipChord = self.tipChord
        sweepAngle = self.sweepAngle
        Xvals = []
        Yvals = []
        Xvals.append(SubComponentPos)
        Yvals.append(BodyTubeDiameter/2)
        Xvals.append(Xvals[0]-span*math.tan(sweepAngle))
        Yvals.append(Yvals[0]+span)
        Xvals.append(Xvals[1]-tipChord)
        Yvals.append(Yvals[1])
        Xvals.append(Xvals[0]-rootChord)
        Yvals.append(Yvals[0])
        Xvals.append(Xvals[0])
        Yvals.append(Yvals[0])
        plt.plot(Xvals,Yvals, color = 'k')

        if self.numFins == 4:
            Yvals = [-1*i for i in Yvals]
            plt.plot(Xvals,Yvals, color = 'k')
            Yvals = [0*i for i in Yvals]
            plt.plot(Xvals,Yvals, color = 'k')

        elif self.numFins == 3:
            Yvals = []
            Yvals.append(-BodyTubeDiameter/2*math.sin(math.radians(30)))
            Yvals.append(-(BodyTubeDiameter/2 + span)*math.sin(math.radians(30)))
            Yvals.append(Yvals[1])
            Yvals.append(Yvals[0])
            Yvals.append(Yvals[0])
            plt.plot (Xvals, Yvals, color = 'k')

Ancestors

Methods

def getAppliedForce(self, rocketState, time, environment, CG)

Expand source code

def getAppliedForce(self, rocketState, time, environment, CG):
    #### If control system exists, use actuator deflections 1:1 to set fin angles ####
    if self.controlSystem != None:
        # Update fin angles
        for i in range(self.numFins):
            self.finList[i].finAngle = self.actuatorList[i].getDeflection(time)

    #### Pre-calculate common properties for all child Fins ####
    precomputedData = self._getPreComputedFinAeroData(rocketState, environment, CG)

    #### Add up forces from all child Fins ####
    aeroForce = ForceMomentSystem(Vector(0,0,0), self.position)
    for fin in self.finList:
        aeroForce += fin.getAppliedForce(rocketState, time, environment, CG, precomputedData)

    # TODO: Correct for sub/transonic rolling moment fin-fin interference from a high number of fins
    
    # Apply fin-body interference factor to total forces in the normal (X/Y) directions and moments
    totalInterferenceFactor = self.bodyOnFinInterferenceFactor * self.finNumberInterferenceFactor
    aeroForce.force.X *= totalInterferenceFactor
    aeroForce.force.Y *= totalInterferenceFactor
    aeroForce.moment.X *= totalInterferenceFactor
    aeroForce.moment.Y *= totalInterferenceFactor

    return aeroForce

def getChord(self, y)

Returns chord at y, where y==0 is at the root chord, and y==span is at the tip chord

Expand source code

def getChord(self, y):
    ''' Returns chord at y, where y==0 is at the root chord, and y==span is at the tip chord '''
    return (self.tipChord - self.rootChord) * (y/self.span) + self.rootChord

def getLogHeader(self)

Expand source code

def getLogHeader(self):
    header = " {}FX(N) {}FY(N) {}FZ(N) {}MX(Nm) {}MY(Nm) {}MZ(Nm)".format(*[self.name]*6)
    
    if self.controlSystem != None:
        for finNumber in range(1, self.numFins+1):
            header += " {}Fin{}Angle(deg)".format(self.name, finNumber)

    return header

def plotShape(self)

Expand source code

def plotShape(self):
    import matplotlib.pyplot as plt

    BodyTubeDiameter = self.bodyRadius*2
    SubComponentPos = self.position.Z
    span = self.span
    rootChord = self.rootChord
    tipChord = self.tipChord
    sweepAngle = self.sweepAngle
    Xvals = []
    Yvals = []
    Xvals.append(SubComponentPos)
    Yvals.append(BodyTubeDiameter/2)
    Xvals.append(Xvals[0]-span*math.tan(sweepAngle))
    Yvals.append(Yvals[0]+span)
    Xvals.append(Xvals[1]-tipChord)
    Yvals.append(Yvals[1])
    Xvals.append(Xvals[0]-rootChord)
    Yvals.append(Yvals[0])
    Xvals.append(Xvals[0])
    Yvals.append(Yvals[0])
    plt.plot(Xvals,Yvals, color = 'k')

    if self.numFins == 4:
        Yvals = [-1*i for i in Yvals]
        plt.plot(Xvals,Yvals, color = 'k')
        Yvals = [0*i for i in Yvals]
        plt.plot(Xvals,Yvals, color = 'k')

    elif self.numFins == 3:
        Yvals = []
        Yvals.append(-BodyTubeDiameter/2*math.sin(math.radians(30)))
        Yvals.append(-(BodyTubeDiameter/2 + span)*math.sin(math.radians(30)))
        Yvals.append(Yvals[1])
        Yvals.append(Yvals[0])
        Yvals.append(Yvals[0])
        plt.plot (Xvals, Yvals, color = 'k')

def precomputeProperties(self)

Expand source code

def precomputeProperties(self):
    self._calculateSweepAngles()

    ### Compute Other Simple Properties ####
    self.planformArea = self.span * (self.tipChord + self.rootChord) / 2
    self.wettedArea = 2*self.planformArea
    self.aspectRatio = (2*self.span)**2/self.planformArea # Aspect ratio of the wing that would be created by reflecting the fin about its root chord (Fleeman/Niskanen)
    self.stallAngle = 15 * (1 + 1.1/self.aspectRatio) # Compute stall angle from Hoerner, 'Fluid Dynamic Lift' for low aspect ratio wings

    # Get the body radius where the fin is mounted
    stageTopZ = self.stage.position.Z
    distanceFromTopOfStageToTopOfThisFinSet = stageTopZ - self.position.Z
    self.bodyRadius = self.stage.getRadius(distanceFromTopOfStageToTopOfThisFinSet)

    self._calculateInterferenceFactors()
    self._calculateMACProperties()
    self._splitFinsIntoSlices()

    self._precomputeSubsonicFinThicknessDrag()
    self._precomputeCPInterpolationPolynomial()

    # Initialize all the child fins in the self.finList list
    self._initChildFins(self.componentDictReader, self.rocket, self.stage)

Inherited members

ActuatedSystem:
- initializeActuators

class FixedForce (componentDictReader, rocket, stage)

Interface definition for rocket components

A Zero-inertia component that applies a constant ForceMomentSystem to the rocket

Expand source code

class FixedForce(RocketComponent):
    def __init__(self, componentDictReader, rocket, stage):
        ''' A Zero-inertia component that applies a constant ForceMomentSystem to the rocket '''
        self.componentDictReader = componentDictReader
        self.rocket = rocket
        self.stage = stage
        self.name = componentDictReader.getDictName()

        # Object is just a force, inertia is zero
        self.inertia = Inertia(Vector(0,0,0), Vector(0,0,0), 0)

        force = componentDictReader.getVector("force")
        forceLocation = componentDictReader.getVector("position")
        moment = componentDictReader.getVector("moment")

        self.force = ForceMomentSystem(force, forceLocation, moment)
        
    def getInertia(self, time, state):
        return self.inertia

    def getAppliedForce(self, rocketState, time, environment, rocketCG):
        return self.force

Ancestors

RocketComponent
abc.ABC

Methods

def getAppliedForce(self, rocketState, time, environment, rocketCG)

Expand source code

def getAppliedForce(self, rocketState, time, environment, rocketCG):
    return self.force

def getInertia(self, time, state)

Expand source code

def getInertia(self, time, state):
    return self.inertia

class FixedMass (componentDictReader, rocket, stage)

Base class for all fixed-mass rocket components Implements functionality to read/store inertia and position info from sim definition file

Expand source code

class FixedMass(RocketComponent):
    '''
        Base class for all fixed-mass rocket components
        Implements functionality to read/store inertia and position info from sim definition file
    '''
    def __init__(self, componentDictReader, rocket, stage):
        self.rocket = rocket
        self.stage = stage
        self.componentDictReader = componentDictReader
        self.name = componentDictReader.getDictName()
            
        mass = componentDictReader.getFloat("mass")

        # Position in simulation definition is relative to stage position
        self.position = componentDictReader.getVector("position") + stage.position # Store position relative to nosecone here
        # CG in simulation definition is relative to component position
        cg = componentDictReader.getVector("cg") + self.position  # Store cg location relative to nosecone here

        try:
            MOI = componentDictReader.getVector("MOI")
        except:
            MOI = Vector(mass*0.01, mass*0.01, mass*0.01) # Avoid having zero moments of inertia

        self.inertia = Inertia(MOI, cg, mass)
        self.zeroForce = ForceMomentSystem(Vector(0,0,0))

    def getInertia(self, time, state):
        return self.inertia

    def getMass(self, time):
        return self.inertia.mass

    def getCG(self, time):
        return self.inertia.CG

    def getAppliedForce(self, rocketState, time, environment, CG):
        return self.zeroForce

Ancestors

RocketComponent
abc.ABC

Subclasses

Methods

def getAppliedForce(self, rocketState, time, environment, CG)

Expand source code

def getAppliedForce(self, rocketState, time, environment, CG):
    return self.zeroForce

def getCG(self, time)

Expand source code

def getCG(self, time):
    return self.inertia.CG

def getInertia(self, time, state)

Expand source code

def getInertia(self, time, state):
    return self.inertia

def getMass(self, time)

Expand source code

def getMass(self, time):
    return self.inertia.mass

class FractionalJetDamping (componentDictReader, rocket, stage)

A component to model Jet damping as per NASA's Two Stage to Orbit verification case

Expand source code

class FractionalJetDamping(RocketComponent):
    ''' A component to model Jet damping as per NASA's Two Stage to Orbit verification case '''

    # Object is just a force, inertia is zero
    inertia = Inertia(Vector(0,0,0), Vector(0,0,0), 0)

    def __init__(self, componentDictReader, rocket, stage):
        self.rocket = rocket
        self.stage = stage
        self.componentDictReader = componentDictReader
        self.name = componentDictReader.getDictName()
        
        self.dampingFraction = componentDictReader.getFloat("fraction")

    def getAppliedForce(self, rocketState, time, environmentalConditions, rocketCG):
        # Only apply damping force if current stage's engine is firing
            # (Other stage's motors will have different exit planes)
        if time > self.stage.motor.ignitionTime and time < self.stage.engineShutOffTime:
            currentRocketInertia = self.rocket.getInertia(time, rocketState)
            
            # Differentiate rate of MOI change
            dt = 0.001
            nextRocketInertia = self.rocket.getInertia(time+dt, rocketState)            
            MOIChangeRate = (currentRocketInertia.MOI.X - nextRocketInertia.MOI.X) / dt

            dampingFactor = MOIChangeRate * self.dampingFraction
            
            angVel = rocketState.angularVelocity
            dampingMoment = Vector(-angVel.X*dampingFactor, -angVel.Y*dampingFactor, 0)

            return ForceMomentSystem(Vector(0,0,0), moment=dampingMoment)
        else:
            return ForceMomentSystem(Vector(0,0,0))

    def getInertia(self, time, state):
        return self.inertia

Ancestors

RocketComponent
abc.ABC

Class variables

var inertia

Methods

def getAppliedForce(self, rocketState, time, environmentalConditions, rocketCG)

Expand source code

def getAppliedForce(self, rocketState, time, environmentalConditions, rocketCG):
    # Only apply damping force if current stage's engine is firing
        # (Other stage's motors will have different exit planes)
    if time > self.stage.motor.ignitionTime and time < self.stage.engineShutOffTime:
        currentRocketInertia = self.rocket.getInertia(time, rocketState)
        
        # Differentiate rate of MOI change
        dt = 0.001
        nextRocketInertia = self.rocket.getInertia(time+dt, rocketState)            
        MOIChangeRate = (currentRocketInertia.MOI.X - nextRocketInertia.MOI.X) / dt

        dampingFactor = MOIChangeRate * self.dampingFraction
        
        angVel = rocketState.angularVelocity
        dampingMoment = Vector(-angVel.X*dampingFactor, -angVel.Y*dampingFactor, 0)

        return ForceMomentSystem(Vector(0,0,0), moment=dampingMoment)
    else:
        return ForceMomentSystem(Vector(0,0,0))

def getInertia(self, time, state)

Expand source code

def getInertia(self, time, state):
    return self.inertia

class NoseCone (componentDictReader, rocket, stage)

Represent a Tangent Ogive Nosecone

Expand source code

class NoseCone(FixedMass, BodyComponent):
    ''' Represent a Tangent Ogive Nosecone '''

    canConnectToComponentAbove = False # Overrides attribute from BodyComponent -> Nosecone must be at top of rocket

    #### Initialization Functions ####
    def __init__(self, componentDictReader, rocket, stage):
        FixedMass.__init__(self, componentDictReader, rocket, stage)
        
        self.aspectRatio = componentDictReader.getFloat("aspectRatio")
        ''' Length over diameter - aka fineness ratio (of exposed surface area) '''

        self.baseDiameter = componentDictReader.getFloat("baseDiameter")
        ''' Diameter at base of nosecone (m) '''

        self.length = self.aspectRatio * self.baseDiameter
        ''' Length from tip to tail, of external surface (m) '''

        self.shape = componentDictReader.getString("shape")
        ''' Description of nosecone shape - something like 'Ogive' or 'Cone' '''
        
        self.surfaceRoughness = componentDictReader.tryGetFloat("surfaceRoughness", defaultValue=self.rocket.surfaceRoughness)
        ''' (m) '''

        self._precomputeGeometry()

    def _precomputeGeometry(self):
        self.length = self.baseDiameter * self.aspectRatio
        self.volume = self._getVolume()
        self.planformArea = self._getPlanformArea()
        self.wettedArea = self._getWettedArea()

        # Using Barrowman's method, the CP Location is a constant
        self.baseArea = self.baseDiameter**2 * math.pi / 4
        self.CPLocation = AeroFunctions.Barrowman_GetCPLocation(self.length, 0, self.baseArea, self.volume)

        # Precompute coefficients for pressure drag interpolation-------------------------------------------        
        #  Niskanen 3.4 and Appendix B
        self.coneHalfAngle = abs(math.atan((self.baseDiameter/2) / self.length))

        self.SubsonicCdPolyCoeffs = computeSubsonicPolyCoeffs(self.coneHalfAngle)
        self.TransonicCdPolyCoeffs = computeTransonicPolyCoeffs(self.coneHalfAngle)

    def _getVolume(self):

        if self.shape == "tangentOgive":
            #https://www.rocketryforum.com/attachments/nosecone_eqn2-pdf.336812/
            baseRadius = self.baseDiameter / 2
            length = self.baseDiameter * self.aspectRatio

            noseconeRadiusOfCurvature = (baseRadius**2 + length**2) / (2*baseRadius)
        
            return math.pi * ( length * noseconeRadiusOfCurvature**2 - length**3 / 3 - (noseconeRadiusOfCurvature - baseRadius) * 
                noseconeRadiusOfCurvature**2 * (math.asin(length / noseconeRadiusOfCurvature)))
        
        elif self.shape == "cone":
            return math.pi * (self.baseDiameter)**2 * self.length / 3

        elif self.shape == "halfSeries":
            #https://www.rocketryforum.com/attachments/nosecone_eqn2-pdf.336812/
            baseRadius = self.baseDiameter / 2
            length = self.baseDiameter * self.aspectRatio

            return math.pi * baseRadius**2 * length / 2

        else:
            raise TypeError("The nosecone shape {} is not recognized".format(self.shape))

    def _getPlanformArea(self):

        if self.shape == "tangentOgive":
            baseRadius = self.baseDiameter / 2
            length = self.baseDiameter * self.aspectRatio

            noseconeRadiusOfCurvature = (baseRadius**2 + length**2) / (2*baseRadius)

            arcSectorAngle = math.asin(length / noseconeRadiusOfCurvature)
          
            arcSectorArea = arcSectorAngle * noseconeRadiusOfCurvature**2 / 2

            return 2 * (arcSectorArea - (length * (noseconeRadiusOfCurvature - baseRadius) / 2))

        elif self.shape == "cone":
            return self.baseDiameter * self.length / 2

        else:
            raise TypeError("The nosecone shape {} is not recognized".format(self.shape))

    def _getWettedArea(self):

        if self.shape == "tangentOgive":
            #TODO Analytical Solution: https://www.rocketryforum.com/attachments/nosecone_eqn2-pdf.336812/
            baseRadius = self.baseDiameter / 2
            length = self.baseDiameter * self.aspectRatio

            noseconeRadiusOfCurvature = (baseRadius**2 + length**2) / (2*baseRadius)
            a = noseconeRadiusOfCurvature

            numSteps = 2000
            stepSize = length / numSteps
            wettedArea = 0

            for x in range(numSteps + 1):
                wettedArea += 2 * math.pi * ((math.sqrt(a**2 - (length - x*stepSize)**2)) + baseRadius - a) * stepSize

            return wettedArea

        elif self.shape == "cone":
            baseRad = self.baseDiameter/2
            return math.pi * baseRad * math.sqrt(baseRad**2 + self.length**2)

        else:
            raise TypeError("The nosecone shape {} is not recognized".format(self.shape))

    #### Operation Functions ####
    def getAppliedForce(self, rocketState, time, environment, CG):
        Aref = self.rocket.Aref
        Mach = AeroParameters.getMachNumber(rocketState, environment)

        # Normal Force --------------------------------------------------------------------------------
        # TODO: Account for rate of pitch/yaw rotation in AOA calculation? Or do separate Pitch/Yaw damping moments?
        AOA = AeroParameters.getTotalAOA(rocketState, environment)
        normalForceCoefficient = AeroFunctions.Barrowman_GetCN(AOA, Aref, 0, self.baseArea)

        # Drag Force ---------------------------------------------------------------------------------------------
        #### Skin Friction ####
        skinFrictionDragCoefficient, rollDampingMoment = AeroFunctions.getCylindricalSkinFrictionDragCoefficientAndRollDampingMoment(rocketState, environment, self.length, Mach,  \
            self.surfaceRoughness, self.wettedArea, Aref, self.rocket.finenessRatio, self.rocket.fullyTurbulentBL)
        
        #### Pressure ####
        pressureDragCoefficient = self._getCd_Pressure(Mach)

        totalDragCoefficient = pressureDragCoefficient + skinFrictionDragCoefficient

        # Damping moments --------------------------------------------------------------------------------------
        # Roll damping due to skin friction
        dampingMoments = rollDampingMoment

        # Combine forces and return total -------------------------------------------------------------------------------------
        return AeroFunctions.forceFromCdCN(rocketState, environment, totalDragCoefficient, normalForceCoefficient, self.CPLocation, Aref, moment=dampingMoments)

    def _getCd_Pressure(self, Mach):
        if Mach < 1:
            # Niskanen pg. 48 eq. 3.87 - Power law interpolation
            pressureDragCoefficient = self.SubsonicCdPolyCoeffs[0] * Mach**self.SubsonicCdPolyCoeffs[1]

        elif Mach > 1 and Mach < 1.3:
            # Interpolate in transonic region - derived from Niskanen Appendix B, Eqns B.3 - B.6
            pressureDragCoefficient = self.TransonicCdPolyCoeffs[0] + self.TransonicCdPolyCoeffs[1]*Mach +  \
                self.TransonicCdPolyCoeffs[2]*Mach**2 + self.TransonicCdPolyCoeffs[3]*Mach**3
            
        else:
            pressureDragCoefficient = getSupersonicPressureDragCoeff_Hoerner(self.coneHalfAngle, Mach)

        return pressureDragCoefficient * self.baseArea/self.rocket.Aref

    def getMaxDiameter(self):
        return self.baseDiameter

    def getRadius(self, distanceFromTip: float) -> float:
        if self.shape == "cone":
            return (self.baseDiameter/2) * (distanceFromTip/self.length)

        elif self.shape == "tangentOgive":
            curveRadius = ((self.baseDiameter/2)**2 + self.length**2) / (self.baseDiameter)
            yComponent = math.sqrt(curveRadius**2 - (self.length - distanceFromTip)**2)
            return yComponent - (curveRadius - (self.baseDiameter/2))

    def plotShape(self) -> None:
        import matplotlib.pyplot as plt
        
        SubComponentPos = self.position.Z
        length = self.baseDiameter * self.aspectRatio
        ConeRad = self.baseDiameter/2
        Rho = (ConeRad**2 + length**2)/(2*ConeRad) # values for tangent ogive equation
        Xvals = list(np.linspace(0,length, num=100)) # create 50 element tuple from 0 to length

        Yvals = [] # create y list

        for x in Xvals:
            # populate y list
            Yvals.append(self.getRadius(x))
        
        Xvals = [-1*i + SubComponentPos for i in Xvals] # these were positive. Need them to be negative to point correct way

        # top line, body tube to tip
        XvalsTop = list(Xvals) 
        YvalsTop = list(Yvals)
        XvalsTop.reverse()
        YvalsTop.reverse()

        # bottom line, tip to body tube
        XvalsBot = list(Xvals)
        YvalsBot = list(Yvals)
        YvalsBot = [-1*i for i in YvalsBot] # so the yvals are negative

        # put it all together now
        Xvals = XvalsTop + XvalsBot
        Yvals = YvalsTop + YvalsBot
        plt.plot(Xvals,Yvals, color = 'k')

Ancestors

Class variables

var canConnectToComponentAbove

Instance variables

var aspectRatio: Length over diameter - aka fineness ratio (of exposed surface area)
var baseDiameter: Diameter at base of nosecone (m)
var length: Length from tip to tail, of external surface (m)
var shape: Description of nosecone shape - something like 'Ogive' or 'Cone'
var surfaceRoughness: (m)

Methods

def getAppliedForce(self, rocketState, time, environment, CG)

Expand source code

def getAppliedForce(self, rocketState, time, environment, CG):
    Aref = self.rocket.Aref
    Mach = AeroParameters.getMachNumber(rocketState, environment)

    # Normal Force --------------------------------------------------------------------------------
    # TODO: Account for rate of pitch/yaw rotation in AOA calculation? Or do separate Pitch/Yaw damping moments?
    AOA = AeroParameters.getTotalAOA(rocketState, environment)
    normalForceCoefficient = AeroFunctions.Barrowman_GetCN(AOA, Aref, 0, self.baseArea)

    # Drag Force ---------------------------------------------------------------------------------------------
    #### Skin Friction ####
    skinFrictionDragCoefficient, rollDampingMoment = AeroFunctions.getCylindricalSkinFrictionDragCoefficientAndRollDampingMoment(rocketState, environment, self.length, Mach,  \
        self.surfaceRoughness, self.wettedArea, Aref, self.rocket.finenessRatio, self.rocket.fullyTurbulentBL)
    
    #### Pressure ####
    pressureDragCoefficient = self._getCd_Pressure(Mach)

    totalDragCoefficient = pressureDragCoefficient + skinFrictionDragCoefficient

    # Damping moments --------------------------------------------------------------------------------------
    # Roll damping due to skin friction
    dampingMoments = rollDampingMoment

    # Combine forces and return total -------------------------------------------------------------------------------------
    return AeroFunctions.forceFromCdCN(rocketState, environment, totalDragCoefficient, normalForceCoefficient, self.CPLocation, Aref, moment=dampingMoments)

def plotShape(self) ‑> None

Expand source code

def plotShape(self) -> None:
    import matplotlib.pyplot as plt
    
    SubComponentPos = self.position.Z
    length = self.baseDiameter * self.aspectRatio
    ConeRad = self.baseDiameter/2
    Rho = (ConeRad**2 + length**2)/(2*ConeRad) # values for tangent ogive equation
    Xvals = list(np.linspace(0,length, num=100)) # create 50 element tuple from 0 to length

    Yvals = [] # create y list

    for x in Xvals:
        # populate y list
        Yvals.append(self.getRadius(x))
    
    Xvals = [-1*i + SubComponentPos for i in Xvals] # these were positive. Need them to be negative to point correct way

    # top line, body tube to tip
    XvalsTop = list(Xvals) 
    YvalsTop = list(Yvals)
    XvalsTop.reverse()
    YvalsTop.reverse()

    # bottom line, tip to body tube
    XvalsBot = list(Xvals)
    YvalsBot = list(Yvals)
    YvalsBot = [-1*i for i in YvalsBot] # so the yvals are negative

    # put it all together now
    Xvals = XvalsTop + XvalsBot
    Yvals = YvalsTop + YvalsBot
    plt.plot(Xvals,Yvals, color = 'k')

Inherited members

BodyComponent:
- getBottomInterfaceLocation
- getMaxDiameter
- getRadius
- getTopInterfaceLocation

class PlanarInterface (location: Vector, component1: RocketComponent, component2: RocketComponent, planeNormal=<MAPLEAF.Motion.CythonVector.Vector object>)

Defines a planar interface between two components
In the local frame, the normalVector is expected to point across the interface from component1 to component2
In a rocket, this means that (with the default normalVector pointing in the -'ve Z direction (towards the tail)), component1 is above component2

Expand source code

class PlanarInterface():
    def __init__(self, location: Vector, component1: RocketComponent, component2: RocketComponent, planeNormal=Vector(0,0,-1)):
        ''' 
            Defines a planar interface between two components  
            In the local frame, the normalVector is expected to point across the interface from component1 to component2   
            In a rocket, this means that (with the default normalVector pointing in the -'ve Z direction (towards the tail)), component1 is above component2   
        '''
        self.location = location
        self.component1 = component1
        self.component2 = component2
        self.normalVector = planeNormal

    @classmethod
    def createPlanarComponentInterfaces(cls, components):
        '''
            Expects components in the list to be sorted by z location (top to bottom)
            Tries to construct PlanarInterface objects connecting all of the BodyComponent object from top to bottom
            Returns a list of PlanarInterface objects, ordered from top to bottom
        '''
        return None
        # Ignore components that aren't of type 'BodyComponent'
        bodyComponents = []
        for comp in components:
            if isinstance(comp, BodyComponent):
                bodyComponents.append(comp)
        
        # Construct interfaces between components
        componentInterfaces = []
        interfaceLocationTolerance = 0.001 # m
        
        for i in range(len(bodyComponents)-1):
            topComponent = bodyComponents[i]
            bottomComponent = bodyComponents[i+1]
            topInterfaceLoc = topComponent.getBottomInterfaceLocation()
            bottomInterfaceLoc = bottomComponent.getTopInterfaceLocation()

            if (topInterfaceLoc - bottomInterfaceLoc).length() < interfaceLocationTolerance:
                interfaceLocation = (topInterfaceLoc + bottomInterfaceLoc) / 2 # Average location is where the interface will be
                componentInterfaces.append(PlanarInterface(interfaceLocation, topComponent, bottomComponent))
            else:
                raise ValueError("Body Component Location mismatch {} ({}) bottom interface at {} vs {} ({}) top interface at {}. Current interface tolerance = 0.001m".format(\
                    topComponent.name, type(topComponent), topInterfaceLoc, bottomComponent.name, type(bottomComponent), bottomInterfaceLoc))

        return componentInterfaces

    @classmethod
    def sortByZLocation(cls, components, state) -> List[RocketComponent]:
        ''' 
            Sort the components in order from top to bottom, component.position.Z
            This function could be relocated somewhere more suitable, at the time of writing, it is only being used to order components before creating interfaces b/w them
        '''
        def getZPosition(component):
            try:
                return component.position.Z
            except AttributeError:
                return component.getInertia(0, state).CG.Z

        components.sort(key=getZPosition, reverse=True)
        return components

Static methods

def createPlanarComponentInterfaces(components)

Expects components in the list to be sorted by z location (top to bottom) Tries to construct PlanarInterface objects connecting all of the BodyComponent object from top to bottom Returns a list of PlanarInterface objects, ordered from top to bottom

Expand source code

@classmethod
def createPlanarComponentInterfaces(cls, components):
    '''
        Expects components in the list to be sorted by z location (top to bottom)
        Tries to construct PlanarInterface objects connecting all of the BodyComponent object from top to bottom
        Returns a list of PlanarInterface objects, ordered from top to bottom
    '''
    return None
    # Ignore components that aren't of type 'BodyComponent'
    bodyComponents = []
    for comp in components:
        if isinstance(comp, BodyComponent):
            bodyComponents.append(comp)
    
    # Construct interfaces between components
    componentInterfaces = []
    interfaceLocationTolerance = 0.001 # m
    
    for i in range(len(bodyComponents)-1):
        topComponent = bodyComponents[i]
        bottomComponent = bodyComponents[i+1]
        topInterfaceLoc = topComponent.getBottomInterfaceLocation()
        bottomInterfaceLoc = bottomComponent.getTopInterfaceLocation()

        if (topInterfaceLoc - bottomInterfaceLoc).length() < interfaceLocationTolerance:
            interfaceLocation = (topInterfaceLoc + bottomInterfaceLoc) / 2 # Average location is where the interface will be
            componentInterfaces.append(PlanarInterface(interfaceLocation, topComponent, bottomComponent))
        else:
            raise ValueError("Body Component Location mismatch {} ({}) bottom interface at {} vs {} ({}) top interface at {}. Current interface tolerance = 0.001m".format(\
                topComponent.name, type(topComponent), topInterfaceLoc, bottomComponent.name, type(bottomComponent), bottomInterfaceLoc))

    return componentInterfaces

def sortByZLocation(components, state) ‑> List[RocketComponent]

Sort the components in order from top to bottom, component.position.Z This function could be relocated somewhere more suitable, at the time of writing, it is only being used to order components before creating interfaces b/w them

Expand source code

@classmethod
def sortByZLocation(cls, components, state) -> List[RocketComponent]:
    ''' 
        Sort the components in order from top to bottom, component.position.Z
        This function could be relocated somewhere more suitable, at the time of writing, it is only being used to order components before creating interfaces b/w them
    '''
    def getZPosition(component):
        try:
            return component.position.Z
        except AttributeError:
            return component.getInertia(0, state).CG.Z

    components.sort(key=getZPosition, reverse=True)
    return components

class RecoverySystem (componentDictReader, rocket, stage)

Represents a recovery system with an arbitrary number of stages

Expand source code

class RecoverySystem(FixedMass):
    ''' Represents a recovery system with an arbitrary number of stages '''

    #### Init Functions ####
    def __init__(self, componentDictReader, rocket, stage):
        FixedMass.__init__(self, componentDictReader, rocket, stage)
        
        self.numStages = componentDictReader.getInt("numStages")

        # Set rocket properties
        rocket.recoverySystem = self

        self._initRecoveryStages()

    def _initRecoveryStages(self):
        self.stageTriggers = [ None, ]
        self.stageTriggerValues = [ None, ]
        self.chuteAreas = [ 0, ]
        self.chuteCds = [ 0, ]
        self.delayTimes = [ 0, ]

        for i in range(1,self.numStages+1):
            self.stageTriggers.append(self.componentDictReader.getString("stage{}Trigger".format(i)))
            try:
                self.stageTriggerValues.append(self.componentDictReader.getFloat("stage{}TriggerValue".format(i)))
            except KeyError:
                self.stageTriggerValues.append(None)
                
            self.chuteAreas.append(self.componentDictReader.getFloat("stage{}ChuteArea".format(i)))
            self.chuteCds.append(self.componentDictReader.getFloat("stage{}Cd".format(i)))
            self.delayTimes.append(self.componentDictReader.getFloat("stage{}DelayTime".format(i)))

        self.currentStage = -1
        self.nextStageDeployTime = None

        # Initialize by deploying stage 0 (no chutes)
        # Sets up trigger function for the real first chute
        self._deployNextStage()

    #### Operational Functions ####
    def getAppliedForce(self, rocketState, time, environment, CG):
        '''
            Calculates force/moment applied by the recovery system using a simple drag coefficient + area model
        '''
        #### Calculate Aero Force ####
        if self.currentStage == 0:
            # No recovery system deployed yet
            return ForceMomentSystem(Vector(0,0,0))
        else:
            # 3DoF force-only aero
            airVel = AeroParameters.getAirVelRelativeToVehicle(rocketState, environment)
            dragForceMagnitude = self.chuteCds[self.currentStage] * self.chuteAreas[self.currentStage] * AeroParameters.getDynamicPressure(rocketState, environment)
            totalForce = airVel.normalize() * dragForceMagnitude
            return ForceMomentSystem(totalForce)

    def _deployNextStage(self):
        ''' 
            Deploys the next stage of the recovery system  
            If appropriate, will set the simulation mode to 3DoF.
            If appropriate, will set up the trigger to deploy the next stage of the recovery system, by calling this function again.
            If this is the final recovery system stage, records the main chute deployment time.
        '''
        self.currentStage += 1
        if self.currentStage > 0:
            print("Deployed Recovery Stage {}".format(self.currentStage))
            
        if self.currentStage == 1:
            self.rocket.isUnderChute = True
            self.rocket._switchTo3DoF()

        if self.currentStage == self.numStages:
            self.rocket.mainChuteDeployTime = self.rocket.rigidBody.time

        if self.currentStage < self.numStages:
            nextTrigger = self.stageTriggers[self.currentStage + 1]
            nextTriggerValue = self.stageTriggerValues[self.currentStage + 1]
            nextTriggerDelay = self.delayTimes[self.currentStage + 1]

            # Translate trigger conditions into language simEventDetector can understand
            if nextTrigger == "Apogee":
                nextTrigger = EventTypes.Apogee
            if nextTrigger == "Altitude":
                nextTrigger = EventTypes.DescendingThroughAltitude
            if nextTrigger == "Time":
                nextTrigger = EventTypes.TimeReached

            # Set the trigger
            self.rocket.simEventDetector.subscribeToEvent(nextTrigger, self._deployNextStage, eventTriggerValue=nextTriggerValue, triggerDelay=nextTriggerDelay)

Ancestors

Methods

def getAppliedForce(self, rocketState, time, environment, CG)

Calculates force/moment applied by the recovery system using a simple drag coefficient + area model

Expand source code

def getAppliedForce(self, rocketState, time, environment, CG):
    '''
        Calculates force/moment applied by the recovery system using a simple drag coefficient + area model
    '''
    #### Calculate Aero Force ####
    if self.currentStage == 0:
        # No recovery system deployed yet
        return ForceMomentSystem(Vector(0,0,0))
    else:
        # 3DoF force-only aero
        airVel = AeroParameters.getAirVelRelativeToVehicle(rocketState, environment)
        dragForceMagnitude = self.chuteCds[self.currentStage] * self.chuteAreas[self.currentStage] * AeroParameters.getDynamicPressure(rocketState, environment)
        totalForce = airVel.normalize() * dragForceMagnitude
        return ForceMomentSystem(totalForce)

class Rocket (rocketDictReader, silent=False, stageToInitialize=None, simRunner=None, environment=None)

Class used to represent a single flying rigid body composed of Stage objects. New instances of this class are also created to model the flight of dropped stages.

Initialization of Rocket(s) is most easily completed through an instance of Simulation To get a single Rocket object, initialize a Simulation and call Simulation.createRocket().
This will return a Rocket initialized on the pad with all its stages, ready for flight.

If initializing manually, can either provide fileName or simDefinition. If a simDefinition is provided, it will be used and fileName will be ignored.

Inputs:

rocketDictReader: (SubDictReader) SubDictReader pointed at the "Rocket" dictionary of the desired simulation definition file.
silent: (bool) controls console output
stageToInitialize: (int or None) controls whether to initialize a complete Rocket or a single (usually dropped) stage. None = initialize complete rocket. n = initialize only stage n, where n >= 1.
simRunner: (Simulation) reference to the current simulation driver/runner
environment: (Environment) environment model from which the rocket will retrieve atmospheric properties and wind speeds

Expand source code

class Rocket(CompositeObject):
    '''
        Class used to represent a single flying rigid body composed of `MAPLEAF.Rocket.Stage` objects.
        New instances of this class are also created to model the flight of dropped stages.
    '''
    #### Initialization ####
    def __init__(self, rocketDictReader, silent=False, stageToInitialize=None, simRunner=None, environment=None):
        '''
            Initialization of Rocket(s) is most easily completed through an instance of Simulation
            To get a single Rocket object, initialize a Simulation and call `MAPLEAF.SimulationRunners.Simulation.createRocket()`.  
            This will return a Rocket initialized on the pad with all its stages, ready for flight.

            If initializing manually, can either provide fileName or simDefinition. If a simDefinition is provided, it will be used and fileName will be ignored.

            Inputs:

            * rocketDictReader:     (`MAPLEAF.IO.SubDictReader`) SubDictReader pointed at the "Rocket" dictionary of the desired simulation definition file.  
            * silent:               (bool) controls console output  
            * stageToInitialize:    (int or None) controls whether to initialize a complete Rocket or a single (usually dropped) stage. None = initialize complete rocket. n = initialize only stage n, where n >= 1.  
            * simRunner:            (`MAPLEAF.SimulationRunners.Simulation`) reference to the current simulation driver/runner
            * environment:          (`MAPLEAF.ENV.Environment`) environment model from which the rocket will retrieve atmospheric properties and wind speeds
        '''
        self.rocketDictReader = rocketDictReader
        self.simDefinition = rocketDictReader.simDefinition

        self.simRunner = simRunner
        ''' Parent instance of `MAPLEAF.SimulationRunners.Simulation` (or derivative sim runner). This is usually the object that has created the current instance of Rocket. '''
        
        self.environment = environment
        ''' Instance of `MAPLEAF.ENV.Environment` '''
        if self.environment == None:
            # If no environment is passed in, create one
            self.environment = Environment(self.simDefinition, silent=silent)

        self.name = rocketDictReader.getString("name")

        self.silent = silent
        ''' Controls output to console '''

        self.stage = stageToInitialize
        ''' If controls whether the whole rocket is initialized (if == None), or a single stage is initialized (Integer stage number) '''

        self.stages = []
        '''
            A list of `MAPLEAF.Rocket.Stage` objects that make up the rocket, ordered from top to bottom.  
            Populated by `_initializeStages`.
        '''

        self.recoverySystem = None
        '''
            Reference to the current Rocket's (which can represent a dropped stage) recovery system. Only a single recovery system is allowed per stage.  
            Set in `MAPLEAF.Rocket.RecoverySystem.__init__`
        '''

        self.rigidBody = None            
        ''' 
            (`MAPLEAF.Motion.RigidBody` or `MAPLEAF.Motion.RigidBody_3DoF`) Responsible for motion integration.  
            Set in `_initializeRigidBody()`.
        '''

        self.isUnderChute = False
        ''' (bool) Controlled by `MAPLEAF.Rocket.Recovery.RecoverySystem._deployNextStage()` '''

        self.mainChuteDeployTime = None
        ''' (float) Filled in during flight by `MAPLEAF.Rocket.Recovery.RecoverySystem._deployNextStage()`  '''
        
        self.engineShutOffTime = None
        ''' (float) Filled in by `MAPLEAF.Rocket.Propulsion.TabulatedMotor.__init__()` upon initialization '''

        self.turbulenceOffWhenUnderChute = rocketDictReader.getBool("Environment.turbulenceOffWhenUnderChute")
        ''' (bool) '''

        self.maxDiameter = self._getMaxBodyTubeDiameter()     
        ''' (float) Holds maximum constant-size body tube diameter, from bodytube components in stages '''

        self.Aref = math.pi * self.maxDiameter**2 / 4
        ''' 
            Reference area for force and moment coefficients.
            Maximum rocket cross-sectional area. Remains constant during flight to retain a 1:1 relationship b/w coefficients in different parts of flight.
            Always based on the maximum body tube diameter in the fully-assembled rocket.
        '''
        
        # TODO: Remove
        self.targetLocation = None

        self.simEventDetector = SimEventDetector(self) 
        ''' (`MAPLEAF.Rocket.SimEventDetector`) Used to trigger things like recovery systems and staging '''

        self.eventTimeStep = rocketDictReader.getFloat("SimControl.TimeStepAdaptation.eventTimingAccuracy")
        ''' If using an adaptive time stepping method, the time step will be overridden near non-time-deterministic discrete events, possibly all the way down to this minimum value '''

        self.addZeroLengthBoatTailsToAccountForBaseDrag = rocketDictReader.getBool("Aero.addZeroLengthBoatTailsToAccountForBaseDrag")
        ''' Controls whether zero-length boat tails are automatically added to the bottom of rockets without them, to make sure base drag is accounted for '''

        self.fullyTurbulentBL = rocketDictReader.getBool("Aero.fullyTurbulentBL")
        ''' Controls whether skin friction is solved assuming a fully turbulent BL or using laminar/transitional flow at lower Reynolds numbers '''

        self.surfaceRoughness = rocketDictReader.getFloat("Aero.surfaceRoughness")
        ''' Default surface roughness for all rocket components '''

        self.finenessRatio = None
        ''' Used in some aerodynamic functions. Updated after initializing subcomponents, and throughout the flight. None if no BodyComponent(s) are present in the rocket '''
        
        self.engineShutOff = False
        '''Used to shut off engines in MAPLEAF.Rocket.Propulsion.DefinedMotor class. Currently set in MAPLE_AF.GNC.Navigation'''

        #### Init Hardware in the loop ####
        subDicts = rocketDictReader.getImmediateSubDicts()
        if "Rocket.HIL" in subDicts:
            self.hardwareInTheLoopControl = "yes"
            quatUpdateRate = rocketDictReader.getInt("HIL.quatUpdateRate")
            posUpdateRate = rocketDictReader.getInt("HIL.posUpdateRate")
            velUpdateRate = rocketDictReader.getInt("HIL.velUpdateRate")
            teensyComPort = rocketDictReader.getString("HIL.teensyComPort")
            imuComPort = rocketDictReader.getString("HIL.imuComPort")
            teensyBaudrate = rocketDictReader.getInt("HIL.teensyBaudrate")
            imuBaudrate = rocketDictReader.getInt("HIL.imuBaudrate")
            self.hilInterface = HILInterface(quatUpdateRate,posUpdateRate,velUpdateRate, teensyComPort, imuComPort, teensyBaudrate, imuBaudrate)
        else:
            self.hardwareInTheLoopControl = "no"

        #### Initialize Logging ####
        self.timeStepLog = None
        ''' Log containing one entry per time step, logs rocket state. None if logging level == 0 '''
        self.derivativeEvaluationLog = None
        ''' Log containing one entry per rocket motion derivative evaluation, contains component forces. None if logging level < 2 '''

        loggingLevel = int(self.simDefinition.getValue("SimControl.loggingLevel"))
        self.loggingLevel = loggingLevel

        if loggingLevel > 0:
            # Create the time step log and add columns to track the rocket state between each time step
            self.timeStepLog = TimeStepLog()
            zeroVector = Vector(0,0,0)
            self.timeStepLog.addColumn("Position(m)", zeroVector)
            self.timeStepLog.addColumn("Velocity(m/s)", zeroVector)
            self.timeStepLog.addColumn("OrientationQuaternion", Quaternion(0,0,0,0))
            self.timeStepLog.addColumn("EulerAngle(rad)", zeroVector)
            self.timeStepLog.addColumn("AngularVelocity(rad/s)", zeroVector)

            if "Adapt" in self.simDefinition.getValue("SimControl.timeDiscretization"):
                self.timeStepLog.addColumn("EstimatedIntegrationError", 0)

        if loggingLevel > 1:
            self.derivativeEvaluationLog = Log()
            zeroVector = Vector(0,0,0)

            self.derivativeEvaluationLog.addColumn("Position(m)", zeroVector)
            self.derivativeEvaluationLog.addColumn("Velocity(m/s)", zeroVector)
            self.derivativeEvaluationLog.addColumn("OrientationQuaternion", Quaternion(0,0,0,0))
            self.derivativeEvaluationLog.addColumn("AngularVelocity(rad/s)", zeroVector)

            self.derivativeEvaluationLog.addColumn("CG(m)", zeroVector)
            self.derivativeEvaluationLog.addColumn("MOI(kg*m^2)", zeroVector)
            self.derivativeEvaluationLog.addColumn("Mass(kg)", 0)
            
            self.derivativeEvaluationLog.addColumn("Wind(m/s)", zeroVector)
            self.derivativeEvaluationLog.addColumn("AirDensity(kg/m^3)", 0)
            self.derivativeEvaluationLog.addColumn("Mach", 0)
            self.derivativeEvaluationLog.addColumn("UnitRe", 0)
            self.derivativeEvaluationLog.addColumn("AOA(deg)", 0)
            self.derivativeEvaluationLog.addColumn("RollAngle(deg)", 0)

            self.derivativeEvaluationLog.addColumn("CPZ(m)", 0)
            self.derivativeEvaluationLog.addColumn("AeroF(N)", zeroVector)
            self.derivativeEvaluationLog.addColumn("AeroM(Nm)", zeroVector)
            self.derivativeEvaluationLog.addColumn("GravityF(N)", zeroVector)
            self.derivativeEvaluationLog.addColumn("TotalF(N)", zeroVector)
        
        #### Init Components ####
        self._createStages()
        self._initializeRigidBody()
        self._switchToStatefulRigidBodyIfRequired()
        self._sortStagesAndComponents()
        self._initializeStagingTriggers()
        self._precomputeComponentProperties()

        #### Init Parent classes ####
        CompositeObject.__init__(self, self.stages)
        self._updateFinenessRatio()

        # Adjust rigid body state position to correspond to the rocket's CG instead of nose cone tip position
        self._moveStatePositionToCG()

        #### Init Guidance/Navigation/Control System (if required) ####
        self.controlSystem = None
        ''' None for uncontrolled rockets. `MAPLEAF.GNC.ControlSystems.RocketControlSystem` for controlled rockets '''
        if ( rocketDictReader.tryGetString("ControlSystem.controlledSystem") != None or rocketDictReader.tryGetString("ControlSystem.MomentController.Type") == "IdealMomentController") and stageToInitialize == None:
            # Only create a control system if this is NOT a dropped stage
            ControlSystemDictReader = SubDictReader("Rocket.ControlSystem", simDefinition=self.simDefinition)
            controlSystemLogging = loggingLevel > 3
            self.controlSystem = RocketControlSystem(ControlSystemDictReader, self, log=controlSystemLogging, silent=silent)

    def _getMaxBodyTubeDiameter(self):
        ''' Gets max body tube diameter directly from config file '''
        stageDicts = self._getStageSubDicts()

        maxDiameter = 0
        for stageDict in stageDicts:
            componentDicts = self.rocketDictReader.getImmediateSubDicts(stageDict)
            
            for componentDict in componentDicts:
                className = self.rocketDictReader.getString(componentDict + ".class")
                
                if className == "Bodytube":
                    diameter = self.rocketDictReader.getFloat(componentDict + ".outerDiameter")
                    maxDiameter = max(maxDiameter, diameter)
        
        return maxDiameter

    def _getStageSubDicts(self):
        # Get all immediate subdictionaries of 'Rocket'
        stageDicts = self.rocketDictReader.getImmediateSubDicts()

        # Assume all subdictionaries represent rocket stages except for these exceptions
        nonStageSubDicts = [ "Rocket.ControlSystem", "Rocket.HIL", "Rocket.Aero" ]
        # Remove them from the list if they're in it
        for dictName in nonStageSubDicts:
            if dictName in stageDicts:
                stageDicts.remove(dictName)

        return stageDicts

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

        initAngVel = AngularVelocity(rotationVector=self.rocketDictReader.getVector("angularVelocity"))
     
        initState_launchTowerFrame = RigidBodyState(initPos, initVel, initOrientation, initAngVel)

        # Convert to the global inertial frame
        initState_globalInertialFrame = self.environment.convertInitialStateToGlobalFrame(initState_launchTowerFrame)

        # Get desired time discretization method
        timeDisc = self.rocketDictReader.getString("SimControl.timeDiscretization")

        #TODO: Check for additional parameters to integrate - if they exist create a StatefulRigidBody + RocketState instead!
            # Ask each rocket component whether it would like to add parameters to integrate after all the components have been initialized

        discardDtCallback = None if (self.simRunner == None) else self.simRunner.discardForceLogsForPreviousTimeStep

        #### Initialize the rigid body ####
        self.rigidBody = RigidBody(
            initState_globalInertialFrame, 
            self._getAppliedForce, 
            self.getInertia, 
            integrationMethod=timeDisc, 
            discardedTimeStepCallback=discardDtCallback,
            simDefinition=self.simDefinition
        )

    def _createStages(self):
        ''' Initialize each of the stages and all of their subcomponents. '''
        stageDicts = self._getStageSubDicts()

        # If we're initializing a dropped stage, figure out which one
        if self.stage != None:
            stageNumbers = []
            stageNumberSet = set()
            for stageDict in stageDicts:
                stageNumber = self.rocketDictReader.getFloat(stageDict + ".stageNumber")
                stageNumbers.append(stageNumber)
                stageNumberSet.add(stageNumber)
            
            if len(stageNumbers) != len(stageNumberSet):
                raise ValueError("For multi-stage rockets, each stage must have a unique stage number. Set the Rocket.StageName.stageNumber value for each stage. 0 for first stage, 1 for second, etc...")
            
            stageNumbers.sort()
            stageToInitialize = stageNumbers[self.stage]

        # Initialize Stage(s)
        initializingAllStages = (self.stage == None)
        for stageDictionary in stageDicts:
            stageDictReader = SubDictReader(stageDictionary, self.simDefinition)
            stageNumber = stageDictReader.getFloat("stageNumber")

            if initializingAllStages or stageNumber == stageToInitialize:
                newStage = Stage(stageDictReader, self)
                self.stages.append(newStage)

                if newStage.name not in self.__dict__: # Avoid clobbering existing info
                    setattr(self, newStage.name, newStage) # Make stage available as rocket.stageName

    def _sortStagesAndComponents(self):
        # Create Planar Interfaces b/w components inside stages
        for stage in self.stages:
            stage.components = PlanarInterface.sortByZLocation(stage.components, self.rigidBody.state)
            stage.componentInterfaces = PlanarInterface.createPlanarComponentInterfaces(stage.components) 

        # Create planar interfaces b/w stages
        self.stages = PlanarInterface.sortByZLocation(self.stages, self.rigidBody.state)
        self.stageInterfaces = PlanarInterface.createPlanarComponentInterfaces(self.stages)

        if self.maxDiameter > 0:
            # Only run this if we're running a real rocket with body tubes
            self._ensureBaseDragIsAccountedFor()

    def _initializeStagingTriggers(self):
        ''' Set up trigger conditions for staging '''       
        # Self.stage is not passed in if the current instance represents a rocket ready to launch - then we have to set up staging events
        if self.stage == None:
            for stageIndex in range(len(self.stages)):
                stage = self.stages[stageIndex]

                if stage.separationConditionType != 'None':
                    self.simEventDetector.subscribeToEvent(
                        stage.separationConditionType, 
                        self._stageSeparation, 
                        stage.separationConditionValue, 
                        triggerDelay=stage.separationDelay
                    )
                
                if stageIndex != len(self.stages)-1:
                    # Set all upper stage motors to ignite a very long time in the future    
                    stage.motor.updateIgnitionTime(1000000000, fakeValue=True)

        else:
            # Otherwise this rocket object is representing a single dropped stage, and no stage separations are necessary
            # Motor is burned out
            self.stages[0].motor.updateIgnitionTime(-1000000000, fakeValue=True)
    
    def _switchToStatefulRigidBodyIfRequired(self):
        '''
            Query all of the rocket components to see if any of them are stateful by attempting to call their getExtraParametersToIntegrate function
            If the rocket contains stateful components, the rocket state is converted to a StateList and all of these state variables requested by components are added to it
            Otherwise, nothing changes and the rocket state remains a RigidBodyState
        '''
        varNames = []
        initVals = []
        derivativeFuncs = []

        # Query all of the components
        for stage in self.stages:
            for component in stage.components:
                try:
                    paramNames, initValues, derivativeFunctions = component.getExtraParametersToIntegrate()
                    varNames += paramNames
                    initVals += initValues
                    derivativeFuncs += derivativeFunctions
                except AttributeError:
                    pass # No extra parameters to integrate

        # If any of the components are stateful, keep track of their states in the main rocket state
        if len(varNames) > 0:
            if len(varNames) != len(initVals) or len(initVals) != len(derivativeFuncs):
                raise ValueError("ERROR: Mismatch in number of extra parameters names ({}), init values({}), and derivative functions({}) to integrate".format(len(varNames), len(initVals), len(derivativeFuncs)))
            
            # Switch to a stateful rigid body
            initState = self.rigidBody.state
            forceFunc = self.rigidBody.forceFunc
            inertiaFunc = self.rigidBody.inertiaFunc
            integrationMethod = self.rocketDictReader.getString("SimControl.timeDiscretization")
            discardDtCallback = None if (self.simRunner == None) else self.simRunner.discardForceLogsForPreviousTimeStep
            self.rigidBody = StatefulRigidBody(
                initState,
                forceFunc,
                inertiaFunc,
                integrationMethod=integrationMethod,
                discardedTimeStepCallback=discardDtCallback,
                simDefinition=self.simDefinition
            )

            # Add the additional state variables
            for i in range(len(varNames)):
                self.rigidBody.addStateVariable(varNames[i], initVals[i], derivativeFuncs[i])

    def _precomputeComponentProperties(self):
        for stage in self.stages:
            for component in stage.components:
                try:
                    component.precomputeProperties()
                except AttributeError:
                    pass

    def _moveStatePositionToCG(self):
        ''' Moves self.rigidBody.state.position to have it represent the rocket's initial CG position, not the initial nose cone position '''
        initInertia = self.getInertia(0, self.rigidBody.state)
        CGPosition_wrtNoseCone_localFrame = initInertia.CG
        CGPosition__wrtNoseCone_globalFrame = self.rigidBody.state.orientation.rotate(CGPosition_wrtNoseCone_localFrame)
        CGPosition_globalFrame = self.rigidBody.state.position + CGPosition__wrtNoseCone_globalFrame
        self.rigidBody.state.position = CGPosition_globalFrame

    def getLength(self):
        totalLength = 0
        for stage in self.stages:
            try:
                totalLength += stage.getLength()
            except TypeError:
                pass # Stage Length was None - no body components in stage

        return totalLength

    def plotShape(self):
        plt.figure()
        plt.gca().set_aspect('equal')

        rocketInertia = self.getInertia(0, self.rigidBody.state)
        TotalCG = rocketInertia.CG.Z
        TotalMass = rocketInertia.mass
        
        TotalCGplt = plt.plot(TotalCG, 0, color='b', marker='d', label='Total CG', linestyle='None')
        CGsubZ = []
        CGsubY = []

        plt.title('Total Rocket CG: %10.4f m \n Total Rocket Mass: %10.4f Kg' % (TotalCG,TotalMass) )

        for stage in self.stages:
            zCGs, yCGs = stage.plotShape()

            # Add subcomponents CGs to arrays
            CGsubZ += zCGs
            CGsubY += yCGs

        SubCGplt = plt.plot(CGsubZ, CGsubY, color='g', marker='.', label='Subcomponent CG', linestyle='None')
        legendHeight = self.maxDiameter
        plt.legend(loc='upper center', bbox_to_anchor = (0.5,-1.05))
        plt.show()

    #### Stage Separation ####
    def _stageSeparation(self):
        print("Stage {} Separation".format(self.simRunner.stagingIndex + 1))
        # Initialize dropped stage as a new rocket
        self.simRunner.createNewDetachedStage()

        # Drop stage from current rocket
        self._dropStage()
        
        # Ignite next motor (set ignition time to the time of stage separation)
        currentTime = self.rigidBody.time
        self.stages[0].motor.updateIgnitionTime(currentTime)

        self._ensureBaseDragIsAccountedFor()
        self._updateFinenessRatio()

    def _dropStage(self, stageIndex=-1):
        droppedStage = self.stages.pop(stageIndex)
        delattr(self, droppedStage.name)
        self.recomputeFixedMassInertia()

        if self.controlSystem != None:
            # Check whether the controlled system was part of the dropped stage
            if self.controlSystem.controlledSystem in droppedStage.components:
                # If so delete the rocket's control system and remove any control system-induced time stepping modifications
                print("Rocket's controlled system was on the dropped stage. Deactivating control system.")
                self.controlSystem.restoreOriginalTimeStepping()
                self.controlSystem = None

    def _ensureBaseDragIsAccountedFor(self):
        ''' If no BoatTail exists at the bottom of the rocket, adds a zero-length boat tail. This is necessary b/c Boat Tail aero-functions are the ones that account for base drag '''
        boatTailComponentAtBottomOfRocket = False
        bottomStage = self.stages[-1]
        for comp in reversed(bottomStage.components):
            if isinstance(comp, BodyComponent):
                if isinstance(comp, BoatTail):
                    boatTailComponentAtBottomOfRocket = True
                break
        
        if not boatTailComponentAtBottomOfRocket and self.addZeroLengthBoatTailsToAccountForBaseDrag:
            if not self.silent:
                print("Adding zero-length BoatTail to the bottom of current bottom stage ({}) to account for base drag".format(bottomStage.name))
            # Create a zero-length, zero-mass boat tail to account for base drag
            zeroInertia = Inertia(Vector(0,0,0), Vector(0,0,0), 0)
            diameter = self.maxDiameter # TODO: Get the actual bottom-body-tube diameter from a future Stage.getRadius function
            length = 0
            position = bottomStage.getBottomInterfaceLocation()
            boatTail = BoatTail(
                diameter, 
                diameter, 
                length, 
                position, 
                zeroInertia, 
                self, 
                bottomStage, 
                "Auto-AddedZeroLengthBoatTail", 
                self.surfaceRoughness
            )
            initializeForceLogging(boatTail, "FakeRocketName.Auto-AddedZeroLengthBoatTail", self)
            bottomStage.components.append(boatTail)     

    def _updateFinenessRatio(self):
        ''' Updates self.finenessRatio based on current BodyComponents in rocket stages '''
        length = self.getLength()
        maxDiameter = max([ stage.getMaxDiameter() for stage in self.stages ])
        if maxDiameter == 0:
            self.finenessRatio = None
        else:
            self.finenessRatio = length/maxDiameter

    #### Component-buildup method for Force ####
    def _getEnvironmentalConditions(self, time, state):
        env = self.environment.getAirProperties(state.position, time)

        # Neglect turbulent component of wind if required
        if self.isUnderChute and self.turbulenceOffWhenUnderChute:
            env = EnvironmentalConditions(
                env.ASLAltitude,
                env.Temp,
                env.Pressure,
                env.Density,
                env.DynamicViscosity,
                env.MeanWind,
                env.MeanWind,
                Vector(0, 0, 0),
            )
        
        return env

    def _getAppliedForce(self, time, state):
        ''' Get the total force currently being experienced by the rocket, used by self.rigidBody to calculate the rocket's acceleration '''
        ### Precomputations and Logging ###
        environment = self._getEnvironmentalConditions(time, state)               # Get and log current air/wind properties
        rocketInertia = self.getInertia(time, state)                                            # Get and log current rocket inertia

        if self.derivativeEvaluationLog is not None:
            self.derivativeEvaluationLog.newLogRow(time)
            self.derivativeEvaluationLog.logValue("Position(m)", state.position)
            self.derivativeEvaluationLog.logValue("Velocity(m/s)", state.velocity)

        ### Component Forces ###
        if not self.isUnderChute:            
            # Precompute and log
            Mach = AeroParameters.getMachNumber(state, environment)
            unitRe = AeroParameters.getReynoldsNumber(state, environment, 1.0)
            AOA = AeroParameters.getTotalAOA(state, environment)
            rollAngle = AeroParameters.getRollAngle(state, environment)

            if self.derivativeEvaluationLog is not None:
                self.derivativeEvaluationLog.logValue("Mach", Mach)
                self.derivativeEvaluationLog.logValue("UnitRe", unitRe)
                self.derivativeEvaluationLog.logValue("AOA(deg)", math.degrees(AOA))
                self.derivativeEvaluationLog.logValue("RollAngle(deg)", rollAngle)
                self.derivativeEvaluationLog.logValue("OrientationQuaternion", state.orientation)
                self.derivativeEvaluationLog.logValue("AngularVelocity(rad/s)", state.angularVelocity)

            # This function will be the inherited function CompositeObject.getAppliedForce
            componentForces = self.getAppliedForce(state, time, environment, rocketInertia.CG) 

        else:
            # When under chute, neglect forces from other components
            componentForces = self.recoverySystem.getAppliedForce(state, time, environment, Vector(0,0,-1))
            
            # Log the recovery system's applied force (Normally handled in CompositeObject.getAppliedForce)
            if hasattr(self.recoverySystem, "forcesLog"):
                self.recoverySystem.forcesLog.append(componentForces.force)
                self.recoverySystem.momentsLog.append(componentForces.moment)
    
        # Move Force-Moment system to rocket CG
        componentForces = componentForces.getAt(rocketInertia.CG)
           
        ### Gravity ###
        gravityForce = self.environment.getGravityForce(rocketInertia, state)        
        totalForce = componentForces + gravityForce

        ### Launch Rail ###
        totalForce = self.environment.applyLaunchTowerForce(state, time, totalForce)        

        if self.derivativeEvaluationLog is not None:
            self.derivativeEvaluationLog.logValue("Wind(m/s)", environment.Wind)
            self.derivativeEvaluationLog.logValue("AirDensity(kg/m^3)", environment.Density)
            self.derivativeEvaluationLog.logValue("CG(m)", rocketInertia.CG)
            self.derivativeEvaluationLog.logValue("MOI(kg*m^2)", rocketInertia.MOI)
            self.derivativeEvaluationLog.logValue("Mass(kg)", rocketInertia.mass)

            CPZ = AeroFunctions._getCPZ(componentForces)            
            self.derivativeEvaluationLog.logValue("CPZ(m)", CPZ)
            self.derivativeEvaluationLog.logValue("AeroF(N)", componentForces.force)
            self.derivativeEvaluationLog.logValue("AeroM(Nm)", componentForces.moment)
            self.derivativeEvaluationLog.logValue("GravityF(N)", gravityForce.force)
            self.derivativeEvaluationLog.logValue("TotalF(N)", totalForce.force)

        return totalForce
    
    #### Driving / Controlling Simulation ####
    def _logState(self):
        '''
            Logs the initial state of the rocket to the time step log
        '''
        state = self.rigidBody.state
        time = self.rigidBody.time

        if self.timeStepLog is not None:
            self.timeStepLog.newLogRow(time)
            self.timeStepLog.logValue("Position(m)", state.position)
            self.timeStepLog.logValue("Velocity(m/s)", state.velocity)

            try: # 6DoF Mode
                self.timeStepLog.logValue("OrientationQuaternion", state.orientation)
                self.timeStepLog.logValue("AngularVelocity(rad/s)", state.angularVelocity)

                # Also log NED Tait-Bryan 3-2-1 z-y-x Euler Angles if in 6DoF mode
                globalOrientation = state.orientation
                orientationOfNEDFrameInGlobalFrame = self.environment.earthModel.getInertialToNEDFrameRotation(*state.position)
                orientationRelativeToNEDFrame = orientationOfNEDFrameInGlobalFrame.conjugate() * globalOrientation
                eulerAngles = orientationRelativeToNEDFrame.toEulerAngles()
                self.timeStepLog.logValue("EulerAngle(rad)", eulerAngles)

            except AttributeError:
                pass # 3DoF mode

        # Print the current time and altitude to the console
        altitude = self.environment.earthModel.getAltitude(*state.position)
        consoleOutput = "{:<8.4f} {:>6.5f}".format(time, altitude)
        print(consoleOutput)

    def _runControlSystem(self):
        '''
            Attempts to run the rocket control system (only runs if it's time to run again, based on its updated rate) (updating target positions for actuators) 
        '''
        if self.controlSystem != None and not self.isUnderChute:
            state = self.rigidBody.state
            time = self.rigidBody.time
            environment = self._getEnvironmentalConditions(time, state)

            ### Run Control Loop ###
            self.controlSystem.runControlLoopIfRequired(time, state, environment)        

    def timeStep(self, dt: float):
        '''
            Tells the simulation to take a time step of size dt.  

            Usually called by functions like `MAPLEAF.SimulationRunners.Simulation.run()`

            Returns:
                * timeStepAdaptationFactor: (float) indicates the factor by which the adaptive time stepping method would like to change the timestep for the next timestep (1.0 for non-adaptive methods)
                * dt:                       (float) actual size of time step taken. Adaptive methods will override the dt asked for if the predicted error for a time step is over their error threshold.
        '''
        # Stop the rocket from sliding off the bottom of the launch rail
        self.rigidBody.state = self.environment.applyLaunchRailMotionConstraints(self.rigidBody.state, self.rigidBody.time)
        
        # Trigger any events that occurred during the last time step
        estimatedTimeToNextEvent, accuratePrediction = self.simEventDetector.triggerEvents()

        # If required, override time step to accurately resolve upcoming discrete events
        if "Adapt" in self.rigidBody.integrate.method and estimatedTimeToNextEvent < dt:
            if accuratePrediction:
                # For time-deterministic events, just set the time step to ever so slightly past the event
                newDt = estimatedTimeToNextEvent + 1e-5
                print("Rocket + SimEventDetector overriding time step from {} to {} to accurately trigger resolve time-deterministic event.".format(dt, newDt))
                dt = newDt
            else:
                # For time-nondeterministic events, slowly approach the event
                newDt = max(estimatedTimeToNextEvent/1.5, self.eventTimeStep)
                estimatedOccurenceTime = self.rigidBody.time + estimatedTimeToNextEvent
                print("Rocket + SimEventDetector overriding time step from {} to {} to accurately resolve upcoming event. Estimated occurence at: {}".format(dt, newDt, estimatedOccurenceTime))
                dt = newDt
                
        self._runControlSystem()
        self._logState()

        # Take timestep
        integrationResult = self.rigidBody.timeStep(dt)

        # If required, log estimated integration error
        if "Adapt" in self.rigidBody.integrate.method and self.timeStepLog is not None:
            self.timeStepLog.logValue("EstimatedIntegrationError", integrationResult.errorMagEstimate)
            
        return integrationResult
    
    def _switchTo3DoF(self):
        ''' Switch to 3DoF simulation after recovery system is deployed '''
        print("Switching to 3DoF simulation")
        new3DoFState = RigidBodyState_3DoF(self.rigidBody.state.position, self.rigidBody.state.velocity)
        
        # Re-read time discretization in case an adaptive method has been selected while using a fixed-update rate control system - in that case, want to switch back to adaptive time stepping for the recovery (uncontrolled) portion of the flight
        originalTimeDiscretization = self.rocketDictReader.getString("SimControl.timeDiscretization")

        if self.simRunner != None:
            self.rigidBody = RigidBody_3DoF(
                new3DoFState, 
                self._getAppliedForce, 
                self.getMass, 
                startTime=self.rigidBody.time, 
                integrationMethod=originalTimeDiscretization, 
                discardedTimeStepCallback=self.simRunner.discardForceLogsForPreviousTimeStep, 
                simDefinition=self.simDefinition
            )
        else:
            self.rigidBody = RigidBody_3DoF(
                new3DoFState, 
                self._getAppliedForce, 
                self.getMass, 
                startTime=self.rigidBody.time, 
                integrationMethod=originalTimeDiscretization, 
                simDefinition=self.simDefinition
            )

    #### After Simulation ####
    def writeLogsToFile(self, directory="."):
        logfilePaths = []

        if self.timeStepLog is not None:
            rocketName = self.components[0].name # Rocket is named after its top stage
            path = os.path.join(directory, "{}_timeStepLog.csv".format(rocketName))
            
            # Time step log
            if self.timeStepLog.writeToCSV(path):
                logfilePaths.append(path)

            # Derivative evaluation log
            if self.derivativeEvaluationLog is not None:
                path = os.path.join(directory, "{}_derivativeEvaluationLog.csv".format(rocketName))  
                
                if self.derivativeEvaluationLog.writeToCSV(path):
                    logfilePaths.append(path)

                    # Calculate aerodynamic coefficients if desired
                    if self.loggingLevel > 2:
                        expandedLogPath = Logging.postProcessForceEvalLog(path, refArea=self.Aref, refLength=self.maxDiameter)
                        logfilePaths.append(expandedLogPath)

            # Control system log
            if self.controlSystem != None:
                if self.controlSystem.log != None:
                    controlSystemLogPath = self.controlSystem.writeLogsToFile(directory)
                    logfilePaths.append(controlSystemLogPath)

        return logfilePaths

Ancestors

CompositeObject

Instance variables

var Aref: Reference area for force and moment coefficients. Maximum rocket cross-sectional area. Remains constant during flight to retain a 1:1 relationship b/w coefficients in different parts of flight. Always based on the maximum body tube diameter in the fully-assembled rocket.
var addZeroLengthBoatTailsToAccountForBaseDrag: Controls whether zero-length boat tails are automatically added to the bottom of rockets without them, to make sure base drag is accounted for
var controlSystem: None for uncontrolled rockets. RocketControlSystem for controlled rockets
var derivativeEvaluationLog: Log containing one entry per rocket motion derivative evaluation, contains component forces. None if logging level < 2
var engineShutOff: Used to shut off engines in MAPLEAF.Rocket.Propulsion.DefinedMotor class. Currently set in MAPLE_AF.GNC.Navigation
var engineShutOffTime: (float) Filled in by TabulatedMotor upon initialization
var environment: Instance of Environment
var eventTimeStep: If using an adaptive time stepping method, the time step will be overridden near non-time-deterministic discrete events, possibly all the way down to this minimum value
var finenessRatio: Used in some aerodynamic functions. Updated after initializing subcomponents, and throughout the flight. None if no BodyComponent(s) are present in the rocket
var fullyTurbulentBL: Controls whether skin friction is solved assuming a fully turbulent BL or using laminar/transitional flow at lower Reynolds numbers
var isUnderChute: (bool) Controlled by RecoverySystem._deployNextStage()
var mainChuteDeployTime: (float) Filled in during flight by RecoverySystem._deployNextStage()
var maxDiameter: (float) Holds maximum constant-size body tube diameter, from bodytube components in stages
var recoverySystem: Reference to the current Rocket's (which can represent a dropped stage) recovery system. Only a single recovery system is allowed per stage.
Set in RecoverySystem
var rigidBody: (RigidBody or RigidBody_3DoF) Responsible for motion integration.
Set in _initializeRigidBody().
var silent: Controls output to console
var simEventDetector: (SimEventDetector) Used to trigger things like recovery systems and staging
var simRunner: Parent instance of Simulation (or derivative sim runner). This is usually the object that has created the current instance of Rocket.
var stage: If controls whether the whole rocket is initialized (if == None), or a single stage is initialized (Integer stage number)
var stages: A list of Stage objects that make up the rocket, ordered from top to bottom.
Populated by _initializeStages.
var surfaceRoughness: Default surface roughness for all rocket components
var timeStepLog: Log containing one entry per time step, logs rocket state. None if logging level == 0
var turbulenceOffWhenUnderChute: (bool)

Methods

def getLength(self)

Expand source code

def getLength(self):
    totalLength = 0
    for stage in self.stages:
        try:
            totalLength += stage.getLength()
        except TypeError:
            pass # Stage Length was None - no body components in stage

    return totalLength

def plotShape(self)

Expand source code

def plotShape(self):
    plt.figure()
    plt.gca().set_aspect('equal')

    rocketInertia = self.getInertia(0, self.rigidBody.state)
    TotalCG = rocketInertia.CG.Z
    TotalMass = rocketInertia.mass
    
    TotalCGplt = plt.plot(TotalCG, 0, color='b', marker='d', label='Total CG', linestyle='None')
    CGsubZ = []
    CGsubY = []

    plt.title('Total Rocket CG: %10.4f m \n Total Rocket Mass: %10.4f Kg' % (TotalCG,TotalMass) )

    for stage in self.stages:
        zCGs, yCGs = stage.plotShape()

        # Add subcomponents CGs to arrays
        CGsubZ += zCGs
        CGsubY += yCGs

    SubCGplt = plt.plot(CGsubZ, CGsubY, color='g', marker='.', label='Subcomponent CG', linestyle='None')
    legendHeight = self.maxDiameter
    plt.legend(loc='upper center', bbox_to_anchor = (0.5,-1.05))
    plt.show()

def timeStep(self, dt: float)

Tells the simulation to take a time step of size dt.

Usually called by functions like Simulation.run()

Returns

timeStepAdaptationFactor: (float) indicates the factor by which the adaptive time stepping method would like to change the timestep for the next timestep (1.0 for non-adaptive methods)
dt: (float) actual size of time step taken. Adaptive methods will override the dt asked for if the predicted error for a time step is over their error threshold.

Expand source code

def timeStep(self, dt: float):
    '''
        Tells the simulation to take a time step of size dt.  

        Usually called by functions like `MAPLEAF.SimulationRunners.Simulation.run()`

        Returns:
            * timeStepAdaptationFactor: (float) indicates the factor by which the adaptive time stepping method would like to change the timestep for the next timestep (1.0 for non-adaptive methods)
            * dt:                       (float) actual size of time step taken. Adaptive methods will override the dt asked for if the predicted error for a time step is over their error threshold.
    '''
    # Stop the rocket from sliding off the bottom of the launch rail
    self.rigidBody.state = self.environment.applyLaunchRailMotionConstraints(self.rigidBody.state, self.rigidBody.time)
    
    # Trigger any events that occurred during the last time step
    estimatedTimeToNextEvent, accuratePrediction = self.simEventDetector.triggerEvents()

    # If required, override time step to accurately resolve upcoming discrete events
    if "Adapt" in self.rigidBody.integrate.method and estimatedTimeToNextEvent < dt:
        if accuratePrediction:
            # For time-deterministic events, just set the time step to ever so slightly past the event
            newDt = estimatedTimeToNextEvent + 1e-5
            print("Rocket + SimEventDetector overriding time step from {} to {} to accurately trigger resolve time-deterministic event.".format(dt, newDt))
            dt = newDt
        else:
            # For time-nondeterministic events, slowly approach the event
            newDt = max(estimatedTimeToNextEvent/1.5, self.eventTimeStep)
            estimatedOccurenceTime = self.rigidBody.time + estimatedTimeToNextEvent
            print("Rocket + SimEventDetector overriding time step from {} to {} to accurately resolve upcoming event. Estimated occurence at: {}".format(dt, newDt, estimatedOccurenceTime))
            dt = newDt
            
    self._runControlSystem()
    self._logState()

    # Take timestep
    integrationResult = self.rigidBody.timeStep(dt)

    # If required, log estimated integration error
    if "Adapt" in self.rigidBody.integrate.method and self.timeStepLog is not None:
        self.timeStepLog.logValue("EstimatedIntegrationError", integrationResult.errorMagEstimate)
        
    return integrationResult

def writeLogsToFile(self, directory='.')

Expand source code

def writeLogsToFile(self, directory="."):
    logfilePaths = []

    if self.timeStepLog is not None:
        rocketName = self.components[0].name # Rocket is named after its top stage
        path = os.path.join(directory, "{}_timeStepLog.csv".format(rocketName))
        
        # Time step log
        if self.timeStepLog.writeToCSV(path):
            logfilePaths.append(path)

        # Derivative evaluation log
        if self.derivativeEvaluationLog is not None:
            path = os.path.join(directory, "{}_derivativeEvaluationLog.csv".format(rocketName))  
            
            if self.derivativeEvaluationLog.writeToCSV(path):
                logfilePaths.append(path)

                # Calculate aerodynamic coefficients if desired
                if self.loggingLevel > 2:
                    expandedLogPath = Logging.postProcessForceEvalLog(path, refArea=self.Aref, refLength=self.maxDiameter)
                    logfilePaths.append(expandedLogPath)

        # Control system log
        if self.controlSystem != None:
            if self.controlSystem.log != None:
                controlSystemLogPath = self.controlSystem.writeLogsToFile(directory)
                logfilePaths.append(controlSystemLogPath)

    return logfilePaths

Inherited members

CompositeObject:
- getMass
- recomputeFixedMassInertia

class RocketComponent (componentDictReader, rocket, stage)

Interface definition for rocket components

Expand source code

class RocketComponent(ABC):
    ''' Interface definition for rocket components '''
    @abstractmethod
    def __init__(self, componentDictReader, rocket, stage):
        return

    @abstractmethod
    def getInertia(self, time, state) -> Inertia:
        return

    @abstractmethod
    def getAppliedForce(self, rocketState, time, environmentalConditions, rocketCG) -> ForceMomentSystem:
        return

    #### Other optional functions ####
    # def getExtraParametersToIntegrate(self):
    '''
        Override this method to make MAPLEAF integrate additional parameters (in addition to rigid body motion)
        Should return three lists
            - A list of (string) parameter names
                - The parameter will then be accessible under state.parameterName whenever the rocket state is available
            - A list of intial parameter values
                - Values can be scalars or vectors, or any other parameter that implements addition, subtraction and multiplication/division by scalars
            - A list of parameter derivative functions
                - Derivative functions should accept two parameters, the current time and the current state object and return the corresponding derivative of the value of that parameter
    '''

Ancestors

abc.ABC

Subclasses

Methods

def getAppliedForce(self, rocketState, time, environmentalConditions, rocketCG) ‑> ForceMomentSystem

Expand source code

@abstractmethod
def getAppliedForce(self, rocketState, time, environmentalConditions, rocketCG) -> ForceMomentSystem:
    return

def getInertia(self, time, state) ‑> Inertia

Expand source code

@abstractmethod
def getInertia(self, time, state) -> Inertia:
    return

class SampleStatefulComponent (componentDictReader, rocket, stage)

Interface definition for rocket components

Expand source code

class SampleStatefulComponent(RocketComponent):
    def __init__(self, componentDictReader, rocket, stage):
        self.rocket = rocket
        self.stage = stage
        self.name = componentDictReader.getDictName()

    def getExtraParametersToIntegrate(self):
        # Examples below for a single parameter to be integrated, can put as many as required in these lists
        paramNames = [ "tankLevel" ]
        initValues = [ 1.0 ]
        derivativeFunctions = [ self.getTankLevelDerivative ]
        
        return paramNames, initValues, derivativeFunctions

    def getTankLevelDerivative(self, time, rocketState):
        return -2*rocketState.tankLevel # tankLevel will asymptotically approach 0

    def getAppliedForce(self, rocketState, time, envConditions, rocketCG):
        mag = -2000*self.getTankLevelDerivative(time, rocketState) # Force magnitude proportional to flow rate out of the tank
        forceVector = Vector(0, 0, mag)

        self.rocket.appendToForceLogLine(" {:>6.4f}".format(mag)) # This will end up in the log file, in the SampleZForce column
        
        return ForceMomentSystem(forceVector)

    def getInertia(self, time, rocketState):
        mass = 5 + rocketState.tankLevel*4.56 # Fixed Mass + fluid mass
        MOI = Vector(mass, mass, mass*0.05) # Related to current mass

        CGz = -3 + rocketState.tankLevel # Moves depending on current tank level
        CG = Vector(0, 0, CGz)
        
        return Inertia(MOI, CG, mass)

Ancestors

RocketComponent
abc.ABC

Methods

def getAppliedForce(self, rocketState, time, envConditions, rocketCG)

Expand source code

def getAppliedForce(self, rocketState, time, envConditions, rocketCG):
    mag = -2000*self.getTankLevelDerivative(time, rocketState) # Force magnitude proportional to flow rate out of the tank
    forceVector = Vector(0, 0, mag)

    self.rocket.appendToForceLogLine(" {:>6.4f}".format(mag)) # This will end up in the log file, in the SampleZForce column
    
    return ForceMomentSystem(forceVector)

def getExtraParametersToIntegrate(self)

Expand source code

def getExtraParametersToIntegrate(self):
    # Examples below for a single parameter to be integrated, can put as many as required in these lists
    paramNames = [ "tankLevel" ]
    initValues = [ 1.0 ]
    derivativeFunctions = [ self.getTankLevelDerivative ]
    
    return paramNames, initValues, derivativeFunctions

def getInertia(self, time, rocketState)

Expand source code

def getInertia(self, time, rocketState):
    mass = 5 + rocketState.tankLevel*4.56 # Fixed Mass + fluid mass
    MOI = Vector(mass, mass, mass*0.05) # Related to current mass

    CGz = -3 + rocketState.tankLevel # Moves depending on current tank level
    CG = Vector(0, 0, CGz)
    
    return Inertia(MOI, CG, mass)

def getTankLevelDerivative(self, time, rocketState)

Expand source code

def getTankLevelDerivative(self, time, rocketState):
    return -2*rocketState.tankLevel # tankLevel will asymptotically approach 0

class SimEventDetector (rocket)

Each sim event detector detects events on a single rocket

Expand source code

class SimEventDetector():
    '''
        Each sim event detector detects events on a single rocket 
    '''

    def __init__(self, rocket):
        self.rocket = rocket

        self.callbackFunctions = []
        self.conditionsEvalFunctions = []
        self.conditionValues = []
        self.triggerDelays = []

    def subscribeToEvent(self, eventType, callbackFunction, eventTriggerValue=None, triggerDelay=0):
        '''
            Inputs:
                eventType:          "apogee", "ascendingThroughAltitude", "descendingThroughAltitude", 'motorBurnout" or "timeReached"
                    Can pass in a member of the EventType Enum defined above, if desired
                callbackfunction:   function to call when the desired event occurs. Must have no arguments
                eventTriggerValue:  None, or a trigger parameter for the eventType - an altitude or a time
                triggerDelay:       (numeric) the amount of time, in seconds, to wait before triggering an event after a trigger condition has been met

            Note:
                Motor Burnout always triggers at the burnout time of the LOWEST stage's motor!
        '''
        # Define map between eventType strings and checker functions
        stringToEvalFnMap = {
            EventTypes.Apogee.value:                    self._isAfterApogee,
            EventTypes.AscendingThroughAltitude.value:  self._isAboveAltitude,
            EventTypes.DescendingThroughAltitude.value: self._isBelowAltitude,
            EventTypes.MotorBurnout.value:              self._isBottomStageMotorBurnedOut,
            EventTypes.TimeReached.value:               self._timeReached
        }

        # Make sure eventType is a string
        if not isinstance(eventType, str):
            eventType = eventType.value # If for example EventType.Apogee is passed in, EventType.Apogee.value retrieves "apogee" from the EventType enum
        
        self.callbackFunctions.append(callbackFunction)
        self.conditionValues.append(eventTriggerValue)
        self.conditionsEvalFunctions.append(stringToEvalFnMap[eventType])
        self.triggerDelays.append(triggerDelay)

    def triggerEvents(self) -> Tuple[float, bool]:
        ''' Checks if any of the events that this SimEventDetector is supposed to detect have happened. If so, triggers their callback functions.
            Returns the estimated time to the next event, and whether that event happens at a set time (time-deterministic) or not (non-time-deterministic - like an altitude condition)
                These return values are used to influence time step adaptation to accurately resolve discrete event timing
        '''
        # Save indices of conditions that have been triggered, will now be removed
        indicesToRemove = []
        
        estimatedTimeToNextEvent = 1e10
        nextEventTimeDeterministic = False
        
        # Precomputed to pass to sub-functions (required for altitude conditions)
        ENUState = self.rocket.environment.convertStateToENUFrame(self.rocket.rigidBody.state)

        # Check each callback condition
        for i in range(len(self.callbackFunctions)):
            # If a conditions has come true
            eventOccurred, timeToOccurrence, timeDeterministic = self.conditionsEvalFunctions[i](self.conditionValues[i], ENUState)

            # Get info about the next event to occur
            if timeToOccurrence < estimatedTimeToNextEvent and not eventOccurred:
                estimatedTimeToNextEvent = timeToOccurrence
                nextEventTimeDeterministic = timeDeterministic

            if eventOccurred:
                if self.triggerDelays[i] == 0:
                    # Call its function
                    # TODO: Print message to simulation log saying that the event has been triggered
                    self.callbackFunctions[i]()
                else:
                    # Schedule the event to happen in triggerDelay seconds
                    triggerTime = self.rocket.rigidBody.time + self.triggerDelays[i]
                    self.subscribeToEvent(EventTypes.TimeReached, self.callbackFunctions[i], triggerTime)

                # And mark it for deletion
                indicesToRemove.append(i)

        # Delete all of the events that were triggered - working in reverse to ensure the indices of the items aren't changed before we delete them
        for i in range(len(indicesToRemove)):
            delI = len(indicesToRemove)-1-i
            indexToRemove = indicesToRemove[delI]
            
            del self.callbackFunctions[indexToRemove]
            del self.conditionValues[indexToRemove]
            del self.conditionsEvalFunctions[indexToRemove]
            del self.triggerDelays[indexToRemove]

        
        # Save velocity and time for next timestep (if we're trying to detect apogee)
            # Doing this here to avoid problems calling the apogee function multiple times in a single time step
        self.lastVelocity = ENUState.velocity
        self.lastTime = self.rocket.rigidBody.time

        return estimatedTimeToNextEvent, nextEventTimeDeterministic
    
    ### Event evaluation functions - each is expected to have a single parameter, return True/False ###
    def _isAfterApogee(self, _, ENUState):
        # Time > 1.0 condition here to avoid setting off events if sliding slightly down before engine lights
        eventOccurred = ENUState.velocity.Z <= 0 and self.rocket.rigidBody.time > 1.0

        timeToOccurrence = 1e10
        try:
            accel = (ENUState.velocity.Z - self.lastVelocity.Z) / (self.rocket.rigidBody.time - self.lastTime)
            if accel < 0 and ENUState.velocity.Z > 0:
                timeToOccurrence = -ENUState.velocity.Z / accel
        except (AttributeError, ZeroDivisionError):
            pass # Haven't saved a velocity/time yet / calling function twice in a single time step

        return eventOccurred, timeToOccurrence, False

    def _isAboveAltitude(self, altitude, ENUState):
        eventOccurred = ENUState.position.Z >= altitude and ENUState.velocity.Z > 0

        if ENUState.velocity.Z > 0:
            timeToOccurrence = (altitude - ENUState.position.Z) / ENUState.velocity.Z
        else:
            timeToOccurrence = 1e10

        return eventOccurred, timeToOccurrence, False

    def _isBelowAltitude(self, altitude, ENUState):
        eventOccurred = ENUState.position.Z <= altitude and ENUState.velocity.Z < 0

        if ENUState.velocity.Z < 0:
            timeToOccurrence = (altitude - ENUState.position.Z) / ENUState.velocity.Z
        else:
            timeToOccurrence = 1e10

        return eventOccurred, timeToOccurrence, False

    def _isBottomStageMotorBurnedOut(self, _, ENUState):
        eventOccurred = self.rocket.rigidBody.time >= self.rocket.stages[-1].engineShutOffTime
        timeToOccurrence = self.rocket.stages[-1].engineShutOffTime - self.rocket.rigidBody.time
        return eventOccurred, timeToOccurrence, True

    def _timeReached(self, time, ENUState):
        eventOccurred = self.rocket.rigidBody.time >= time
        timeToOccurrence = time - self.rocket.rigidBody.time
        return eventOccurred, timeToOccurrence, True

Methods

def subscribeToEvent(self, eventType, callbackFunction, eventTriggerValue=None, triggerDelay=0)

Inputs

eventType: "apogee", "ascendingThroughAltitude", "descendingThroughAltitude", 'motorBurnout" or "timeReached" Can pass in a member of the EventType Enum defined above, if desired callbackfunction: function to call when the desired event occurs. Must have no arguments eventTriggerValue: None, or a trigger parameter for the eventType - an altitude or a time triggerDelay: (numeric) the amount of time, in seconds, to wait before triggering an event after a trigger condition has been met

Note

Motor Burnout always triggers at the burnout time of the LOWEST stage's motor!

Expand source code

def subscribeToEvent(self, eventType, callbackFunction, eventTriggerValue=None, triggerDelay=0):
    '''
        Inputs:
            eventType:          "apogee", "ascendingThroughAltitude", "descendingThroughAltitude", 'motorBurnout" or "timeReached"
                Can pass in a member of the EventType Enum defined above, if desired
            callbackfunction:   function to call when the desired event occurs. Must have no arguments
            eventTriggerValue:  None, or a trigger parameter for the eventType - an altitude or a time
            triggerDelay:       (numeric) the amount of time, in seconds, to wait before triggering an event after a trigger condition has been met

        Note:
            Motor Burnout always triggers at the burnout time of the LOWEST stage's motor!
    '''
    # Define map between eventType strings and checker functions
    stringToEvalFnMap = {
        EventTypes.Apogee.value:                    self._isAfterApogee,
        EventTypes.AscendingThroughAltitude.value:  self._isAboveAltitude,
        EventTypes.DescendingThroughAltitude.value: self._isBelowAltitude,
        EventTypes.MotorBurnout.value:              self._isBottomStageMotorBurnedOut,
        EventTypes.TimeReached.value:               self._timeReached
    }

    # Make sure eventType is a string
    if not isinstance(eventType, str):
        eventType = eventType.value # If for example EventType.Apogee is passed in, EventType.Apogee.value retrieves "apogee" from the EventType enum
    
    self.callbackFunctions.append(callbackFunction)
    self.conditionValues.append(eventTriggerValue)
    self.conditionsEvalFunctions.append(stringToEvalFnMap[eventType])
    self.triggerDelays.append(triggerDelay)

def triggerEvents(self) ‑> Tuple[float, bool]

Checks if any of the events that this SimEventDetector is supposed to detect have happened. If so, triggers their callback functions. Returns the estimated time to the next event, and whether that event happens at a set time (time-deterministic) or not (non-time-deterministic - like an altitude condition) These return values are used to influence time step adaptation to accurately resolve discrete event timing

Expand source code

def triggerEvents(self) -> Tuple[float, bool]:
    ''' Checks if any of the events that this SimEventDetector is supposed to detect have happened. If so, triggers their callback functions.
        Returns the estimated time to the next event, and whether that event happens at a set time (time-deterministic) or not (non-time-deterministic - like an altitude condition)
            These return values are used to influence time step adaptation to accurately resolve discrete event timing
    '''
    # Save indices of conditions that have been triggered, will now be removed
    indicesToRemove = []
    
    estimatedTimeToNextEvent = 1e10
    nextEventTimeDeterministic = False
    
    # Precomputed to pass to sub-functions (required for altitude conditions)
    ENUState = self.rocket.environment.convertStateToENUFrame(self.rocket.rigidBody.state)

    # Check each callback condition
    for i in range(len(self.callbackFunctions)):
        # If a conditions has come true
        eventOccurred, timeToOccurrence, timeDeterministic = self.conditionsEvalFunctions[i](self.conditionValues[i], ENUState)

        # Get info about the next event to occur
        if timeToOccurrence < estimatedTimeToNextEvent and not eventOccurred:
            estimatedTimeToNextEvent = timeToOccurrence
            nextEventTimeDeterministic = timeDeterministic

        if eventOccurred:
            if self.triggerDelays[i] == 0:
                # Call its function
                # TODO: Print message to simulation log saying that the event has been triggered
                self.callbackFunctions[i]()
            else:
                # Schedule the event to happen in triggerDelay seconds
                triggerTime = self.rocket.rigidBody.time + self.triggerDelays[i]
                self.subscribeToEvent(EventTypes.TimeReached, self.callbackFunctions[i], triggerTime)

            # And mark it for deletion
            indicesToRemove.append(i)

    # Delete all of the events that were triggered - working in reverse to ensure the indices of the items aren't changed before we delete them
    for i in range(len(indicesToRemove)):
        delI = len(indicesToRemove)-1-i
        indexToRemove = indicesToRemove[delI]
        
        del self.callbackFunctions[indexToRemove]
        del self.conditionValues[indexToRemove]
        del self.conditionsEvalFunctions[indexToRemove]
        del self.triggerDelays[indexToRemove]

    
    # Save velocity and time for next timestep (if we're trying to detect apogee)
        # Doing this here to avoid problems calling the apogee function multiple times in a single time step
    self.lastVelocity = ENUState.velocity
    self.lastTime = self.rocket.rigidBody.time

    return estimatedTimeToNextEvent, nextEventTimeDeterministic

class Stage (stageDictReader, rocket)

Represents a single rocket stage

Expand source code

class Stage(CompositeObject, BodyComponent):
    ''' Represents a single rocket stage '''

    def __init__(self, stageDictReader, rocket):
        #TODO: This should take a boolean parameter 'dropped' which specified whether the stage being initialized is on the pad (full engine) or being dropped (empty engine)
        self.rocket = rocket
        self.stageDictReader = stageDictReader

        self.name = stageDictReader.getDictName()
        self.stageNumber = stageDictReader.getInt("stageNumber")
        self.position = stageDictReader.getVector("position")

        self.separationConditionType = stageDictReader.getString("separationTriggerType")
        self.separationDelay = stageDictReader.getFloat("separationDelay")
        self.separationConditionValue = stageDictReader.getFloat("separationTriggerValue")

        #### Placeholder fields, to be filled by subcomponents ####
        self.motor = None # Filled in by Motor.__init__
        self.engineShutOffTime = None # Filled in by Motor.__init__ - used for base drag calculations

        self.components = []
        self.initializeSubComponents()

        # Check if any stage-level inertia properties are being overriden
        CGOverride = stageDictReader.tryGetVector("constCG", defaultValue=None)
        if CGOverride != None:
            CGOverride += self.position
        massOverride = stageDictReader.tryGetFloat("constMass", defaultValue=None)
        MOIOverride = stageDictReader.tryGetVector("constMOI", defaultValue=None)

        # Initialize functionality inherited from CompositeObject - pass in any overrides
        CompositeObject.__init__(self, components=self.components, CGOverride=CGOverride, MOIOverride=MOIOverride, massOverride=massOverride)
        
    def initializeSubComponents(self):
        # Assume each sub dictionary represents a component
        subDicts = self.stageDictReader.getImmediateSubDicts()

        def returnZeroIfBodyComponent(componentDictPath):
            className = self.stageDictReader.getString(componentDictPath + ".class")
            componentClass = stringNameToClassMap[className]
            
            if issubclass(componentClass, BodyComponent):
                return 0, componentDictPath
            else: 
                return 1, componentDictPath

        # Make sure BodyComponents are initializede before others
            # This allows non-body components (which are attached to body components) like the fins
                # to get info about the body radius where they're mounted
        subDicts.sort(key=returnZeroIfBodyComponent)

        # Initialize each component, add to componets
        for componentDict in subDicts:
            newComponent = rocketComponentFactory(componentDict, self.rocket, self)
            self.components.append(newComponent)

            try:
                if newComponent.name not in self.__dict__: # Avoid clobbering existing info
                    setattr(self, newComponent.name, newComponent) # Make component available as stage.componentName
            except AttributeError:
                pass # Expect to arrive here if the component doesn't define a "name" attribute   

    def getComponentsOfType(self, type):
        matchingComponents = []
        for comp in self.components:
            if isinstance(comp, type):
                matchingComponents.append(comp)

        return matchingComponents

    #### Geometry / Introspection ####
    def _getFirstComponentOfType(self, componentList, desiredType):
        for comp in componentList:
            if isinstance(comp, desiredType):
                return comp
        return None

    def getTopInterfaceLocation(self) -> Union[Vector, None]:
        # Expects components to be sorted by z-location, top to bottom
        comp = self._getFirstComponentOfType(self.components, BodyComponent)
        try:
            topZ = comp.getTopInterfaceLocation().Z
            return Vector(self.position.X, self.position.Y, topZ)   
        except AttributeError:
            return None # Occurs when comp.getTopInterfaceLocation() == None

    def getBottomInterfaceLocation(self) -> Union[Vector, None]:
        # Expects components to be sorted by z-location, top to bottom
        comp = self._getFirstComponentOfType(reversed(self.components), BodyComponent)
        try:
            bottomZ = comp.getBottomInterfaceLocation().Z
            return Vector(self.position.X, self.position.Y, bottomZ)   
        except AttributeError:
            return None # Occurs when comp.getBottomInterfaceLocation() == None

    def getLength(self):
        try:
            topBodyComponent = self._getFirstComponentOfType(self.components, BodyComponent)
            topLocationZ = topBodyComponent.position.Z

            bottomBodyComponent = self._getFirstComponentOfType(reversed(self.components), BodyComponent)
            bottomLocationZ = bottomBodyComponent.position.Z - bottomBodyComponent.length
            return (topLocationZ - bottomLocationZ)
        except AttributeError:
            return None # No body components in current Stage

    def getMaxDiameter(self):
        maxDiameter = 0
        for comp in self.components:
            try:
                maxDiameter = max(comp.getMaxDiameter(), maxDiameter)
            except AttributeError:
                pass # Component is not a BodyComponent
        
        return maxDiameter

    def getRadius(self, distanceFromTop):
        topOfStageZ = self.position.Z
        desiredZ = topOfStageZ - distanceFromTop

        for comp in self.components:
            try:
                # Check if desired location is part of this component, if so, return the component's radius
                compTopZ = comp.position.Z
                compBottomZ = compTopZ - comp.length
                if desiredZ <= compTopZ and desiredZ >= compBottomZ:
                    distanceFromTopOfComponent = compTopZ - desiredZ
                    return comp.getRadius(distanceFromTopOfComponent)

            except AttributeError:
                pass # Not a body component
        
        raise ValueError("Length {} from the top of stage '{}' does not fall between the top and bottom of any of its components".format(distanceFromTop, self.name))

    def plotShape(self):
        '''
            Calls .plotShape() on all subcomponents
            Returns CGsubZ, CGsubY (two lists of scalars for a 2D plot)
        '''
        CGsubZ = []
        CGsubY = []
        for component in self.components:               
            try:
                component.plotShape()
                initialComponentInertia = component.getInertia(0, self.rocket.rigidBody.state)
                CGsubZ.append(initialComponentInertia.CG.Z)
                CGsubY.append(initialComponentInertia.CG.Y)
            
            except AttributeError:
                pass # Plotting function not implemented
        
        return CGsubZ, CGsubY

Ancestors

Methods

def getComponentsOfType(self, type)

Expand source code

def getComponentsOfType(self, type):
    matchingComponents = []
    for comp in self.components:
        if isinstance(comp, type):
            matchingComponents.append(comp)

    return matchingComponents

def getLength(self)

Expand source code

def getLength(self):
    try:
        topBodyComponent = self._getFirstComponentOfType(self.components, BodyComponent)
        topLocationZ = topBodyComponent.position.Z

        bottomBodyComponent = self._getFirstComponentOfType(reversed(self.components), BodyComponent)
        bottomLocationZ = bottomBodyComponent.position.Z - bottomBodyComponent.length
        return (topLocationZ - bottomLocationZ)
    except AttributeError:
        return None # No body components in current Stage

def initializeSubComponents(self)

Expand source code

def initializeSubComponents(self):
    # Assume each sub dictionary represents a component
    subDicts = self.stageDictReader.getImmediateSubDicts()

    def returnZeroIfBodyComponent(componentDictPath):
        className = self.stageDictReader.getString(componentDictPath + ".class")
        componentClass = stringNameToClassMap[className]
        
        if issubclass(componentClass, BodyComponent):
            return 0, componentDictPath
        else: 
            return 1, componentDictPath

    # Make sure BodyComponents are initializede before others
        # This allows non-body components (which are attached to body components) like the fins
            # to get info about the body radius where they're mounted
    subDicts.sort(key=returnZeroIfBodyComponent)

    # Initialize each component, add to componets
    for componentDict in subDicts:
        newComponent = rocketComponentFactory(componentDict, self.rocket, self)
        self.components.append(newComponent)

        try:
            if newComponent.name not in self.__dict__: # Avoid clobbering existing info
                setattr(self, newComponent.name, newComponent) # Make component available as stage.componentName
        except AttributeError:
            pass # Expect to arrive here if the component doesn't define a "name" attribute

def plotShape(self)

Calls .plotShape() on all subcomponents Returns CGsubZ, CGsubY (two lists of scalars for a 2D plot)

Expand source code

def plotShape(self):
    '''
        Calls .plotShape() on all subcomponents
        Returns CGsubZ, CGsubY (two lists of scalars for a 2D plot)
    '''
    CGsubZ = []
    CGsubY = []
    for component in self.components:               
        try:
            component.plotShape()
            initialComponentInertia = component.getInertia(0, self.rocket.rigidBody.state)
            CGsubZ.append(initialComponentInertia.CG.Z)
            CGsubY.append(initialComponentInertia.CG.Y)
        
        except AttributeError:
            pass # Plotting function not implemented
    
    return CGsubZ, CGsubY

Inherited members

CompositeObject:
- getMass
- recomputeFixedMassInertia
BodyComponent:
- getBottomInterfaceLocation
- getMaxDiameter
- getRadius
- getTopInterfaceLocation

class TabulatedAeroForce (componentDictReader, rocket, stage)

A zero-inertia component with aerodynamic coefficients that are tabulated according to one or more parameters (ex. AOA)

Expand source code

class TabulatedAeroForce(AeroForce):
    ''' A zero-inertia component with aerodynamic coefficients that are tabulated according to one or more parameters (ex. AOA) '''

    def __init__(self, componentDictReader, rocket, stage):
        self.componentDictReader = componentDictReader
        self.rocket = rocket
        self.stage = stage
        self.name = componentDictReader.getDictName()

        self.position = componentDictReader.getVector("position")
        self.Aref = componentDictReader.getFloat("Aref")
        self.Lref = componentDictReader.getFloat("Lref")

        coefficientTableFilePath = componentDictReader.getString("filePath")
        self._loadCoefficients(coefficientTableFilePath)

    def _loadCoefficients(self, filePath):
        # Load first row to figure out what the columns mean
        with open(filePath) as f:
            columnNames = f.readline().strip().split(',')

        # Get functions that calculate the parameters used for interpolation
        # All these 'key'/parameter columns are expected to come before 'value' columns to be interpolated over
        self.parameterFunctions = []
        i = 0
        while i < len(columnNames):
            col = columnNames[i]
            if col in AeroParameters.stringToAeroFunctionMap:
                self.parameterFunctions.append(AeroParameters.stringToAeroFunctionMap[col])
            else:
                break
            i += 1

        # Continue parsing column names - aero coefficient names now            
        # This is the ordering expected by AeroFunctions.forceFromCoefficients
        aeroCoeffStrings = [ "CD", "CL", "CMx", "CMy", "CMz" ]
        self.aeroCoeffIndices = [] # Provides mapping between value column position in interpolation table & position in output aero coefficient list (ordered like aeroCoeffStrings above)
        while i < len(columnNames):
            coeff = columnNames[i]

            if coeff in aeroCoeffStrings:
                self.aeroCoeffIndices.append(aeroCoeffStrings.index(coeff))
                
            else:
                raise ValueError("ERROR: One of the following columns: {} did not match any of the expected columns names: Keys: {}, values: {}. \
                    Or was in the wrong order. All key columns must come BEFORE value columns.".format(columnNames, AeroParameters.stringToAeroFunctionMap.keys(), aeroCoeffStrings))
            i += 1

        # Load the data table to be interpolated
        dataTable = np.loadtxt(filePath, delimiter=',', skiprows=1)

        nKeyCols = len(self.parameterFunctions)
        keys = dataTable[:, 0:nKeyCols]
        aeroCoefficients = dataTable[:, nKeyCols:]

        if nKeyCols > 1:
            # Create n-dimensional interpolation function for aero coefficients
            self._interpAeroCoefficients =  NoNaNLinearNDInterpolator(keys, aeroCoefficients, filePath)
        else:
            # Save to use with MAPLEAF.Motion.linInterp
            self.keys = [ key[0] for key in keys ]
            self.values = aeroCoefficients

    def _getAeroCoefficients(self, state, environment):
        keys = AeroParameters.getAeroPropertiesList(self.parameterFunctions, state, environment)

        if len(keys) > 1:
            # Multi-dimensional linear interpolation
            interpolatedCoefficients = self._interpAeroCoefficients(keys)[0]
        else:
            # 1D linear interpolation
            interpolatedCoefficients = linInterp(self.keys, self.values, keys[0])

        aeroCoefficients = [0.0] * 5
        for i in range(len(interpolatedCoefficients)):
            indexInCoeffArray = self.aeroCoeffIndices[i]
            aeroCoefficients[indexInCoeffArray] = interpolatedCoefficients[i]

        return aeroCoefficients

    def getAppliedForce(self, state, time, environment, rocketCG):
        aeroCoefficients = self._getAeroCoefficients(state, environment)
        return AeroFunctions.forceFromCoefficients(state, environment, *aeroCoefficients, self.position, self.Aref, self.Lref)

Ancestors

Methods

def getAppliedForce(self, state, time, environment, rocketCG)

Expand source code

def getAppliedForce(self, state, time, environment, rocketCG):
    aeroCoefficients = self._getAeroCoefficients(state, environment)
    return AeroFunctions.forceFromCoefficients(state, environment, *aeroCoefficients, self.position, self.Aref, self.Lref)

class TabulatedInertia (componentDictReader, rocket, stage)

A zero-force component with time-varying tabulated inertia

Expand source code

class TabulatedInertia(RocketComponent):
    ''' A zero-force component with time-varying tabulated inertia '''
    def __init__(self, componentDictReader, rocket, stage):
        self.rocket = rocket
        self.stage = stage
        self.componentDictReader = componentDictReader
        self.name = componentDictReader.getDictName()

        self.zeroForce = ForceMomentSystem(Vector(0,0,0))

        inertiaTableFilePath = componentDictReader.getString("filePath")
        self._parseInertiaTable(inertiaTableFilePath)

    def _parseInertiaTable(self, filePath):
        data = np.loadtxt(filePath, skiprows=1, delimiter=',')
        self.times = data[:, 0]
        self.inertiaData = data[:, 1:]

        # Check that the right number of columns is present
        if data.shape[1] != 8:
            raise ValueError("Wrong number of columns in inertia table: {}. Expecting 8 columns: \
                Time, Mass, CGx, CGy, CGz, MOIx, MOIy, MOIz")

    def getInertia(self, time, state):
        inertiaData = linInterp(self.times, self.inertiaData, time)
        # MOI is last three columns, CG is the three before that, and mass is column 0
        return Inertia(Vector(*inertiaData[-3:]), Vector(*inertiaData[1:4]), inertiaData[0])
    
    def getAppliedForce(self, rocketState, time, environment, CG):
        return self.zeroForce

Ancestors

RocketComponent
abc.ABC

Methods

def getAppliedForce(self, rocketState, time, environment, CG)

Expand source code

def getAppliedForce(self, rocketState, time, environment, CG):
    return self.zeroForce

def getInertia(self, time, state)

Expand source code

def getInertia(self, time, state):
    inertiaData = linInterp(self.times, self.inertiaData, time)
    # MOI is last three columns, CG is the three before that, and mass is column 0
    return Inertia(Vector(*inertiaData[-3:]), Vector(*inertiaData[1:4]), inertiaData[0])

class TabulatedMotor (componentDictReader, rocket, stage)

Interface

Initialization: In rocket text file, attribute: "path", pointing at a motor definition text file Format:"test/testMotorDefintion.txt" Functions: .Thrust(time) returns current thrust level .OxWeight(time) returns current oxidizer weight .FuelWeight(time) returns current fuel weight Attributes: .initialOxidizerWeight .initialFuelWeight

All in same units as in the motor definition file (linearly interpolated) The motor is assumed to apply zero moment to the rocket, thrust is applied in z-direction

Expand source code

class TabulatedMotor(RocketComponent):
    '''
    Interface:
        Initialization:
            In rocket text file, attribute: "path", pointing at a motor definition text file
            Format:"test/testMotorDefintion.txt"
        Functions:
            .Thrust(time) returns current thrust level
            .OxWeight(time) returns current oxidizer weight
            .FuelWeight(time) returns current fuel weight
        Attributes:
            .initialOxidizerWeight
            .initialFuelWeight

        All in same units as in the motor definition file (linearly interpolated)
        The motor is assumed to apply zero moment to the rocket, thrust is applied in z-direction
    '''

    #### Init Functions ####
    def __init__(self, componentDictReader, rocket, stage):
        #TODO: Oxidizer and Fuel CG Locations should be defined relative to the motor location
        self.rocket = rocket
        self.stage = stage
        self.componentDictReader = componentDictReader
        self.name = componentDictReader.getDictName()

        stage.motor = self
        self.classType = componentDictReader.getString("class")
        self.ignitionTime = 0 # Modified by Rocket._initializeStaging and Rocket._stageSeparation

        # Impulse adjustments (mostly for Monte Carlo sims)
        self.impulseAdjustFactor = componentDictReader.getFloat("impulseAdjustFactor")
        self.burnTimeAdjustFactor = componentDictReader.getFloat("burnTimeAdjustFactor")

        motorFilePath = componentDictReader.getString("path")
        self._parseMotorDefinitionFile(motorFilePath)

        # Set the position to the initial CG location
        initInertia = self.getInertia(0, "fakeState")
        self.position = initInertia.CG

    #TODO: Build converter/parser for standard engine format like rasp/.eng or something like that

    def _parseMotorDefinitionFile(self, motorFilePath):
        ''' Parses a motor definition text file. See MAPLEAF/Examples/Motors for examples '''
      
        # Get motor definition text
        with open(motorFilePath, "r") as motorFile:
            motorFileText = motorFile.read()

        # Remove all comment rows
        comment = re.compile("#.*") 
        motorFileText = re.sub(comment, "", motorFileText)
        
        #Remove blank lines
        motorFileText = [line for line in motorFileText.split('\n') if line.strip() != '']
        
        # Parse CG locations
        # TODO: Future motors should be able to exist off the rocket's center axis
        self.initOxCG_Z = float(motorFileText[0].split()[1]) + self.stage.position.Z
        self.finalOxCG_Z = float(motorFileText[1].split()[1]) + self.stage.position.Z
        self.initFuelCG_Z = float(motorFileText[2].split()[1]) + self.stage.position.Z
        self.finalFuelCG_Z = float(motorFileText[3].split()[1]) + self.stage.position.Z
        motorFileText = motorFileText[4:]

        # Parse data; Columns defined in MAPLEAF/Examples/Motors/test.txt
        # Gets defined values for: Time, thrust, oxFlowRate, fuelFlowRate, oxMOI, fuelMOI
        self.times = []
        self.thrustLevels = []
        oxFlowRate = []
        fuelFlowRate = []
        self.oxMOIs = []
        self.fuelMOIs = []
        for dataLine in motorFileText:
            # Splits line at each white space
            info = dataLine.split()

            self.times.append(float(info[0]))
            self.thrustLevels.append(float(info[1]))
            oxFlowRate.append(float(info[2]))
            fuelFlowRate.append(float(info[3]))
            
            oxVecStartIndex = dataLine.index('(')
            oxVecEndIndex = dataLine.index(')', oxVecStartIndex)+1
            oxVecString  = dataLine[oxVecStartIndex:oxVecEndIndex]
            oxMOIVec = Vector(oxVecString)
            self.oxMOIs.append(oxMOIVec)

            fuelVecStartIndex = dataLine.index('(', oxVecEndIndex)
            fuelVecEndIndex = dataLine.index(')', fuelVecStartIndex)+1
            fuelVecString  = dataLine[fuelVecStartIndex:fuelVecEndIndex]
            fuelMOIVec = Vector(fuelVecString)
            self.fuelMOIs.append(fuelMOIVec)

        # Tell the rocket and stage when their engines shut off -> used for flight animations
        self.stage.engineShutOffTime = self.times[-1]
        if self.rocket.engineShutOffTime == None:
            self.rocket.engineShutOffTime = self.times[-1]
        else:
            self.rocket.engineShutOffTime = max(self.rocket.engineShutOffTime, self.times[-1])

        # Apply adjustment factors for monte carlo sims
        self.thrustLevels = [ thrust*self.impulseAdjustFactor/self.burnTimeAdjustFactor for thrust in self.thrustLevels ]
        self.times = [ t*self.burnTimeAdjustFactor for t in self.times ]

        # Calculate initial fuel and oxidizer masses through trapezoid rule
        # Trapezoid rule matches the linear interpolation used to find thrust values
        self.initialOxidizerWeight = 0
        self.initialFuelWeight = 0
        self.oxWeights = [ 0 ]
        self.fuelWeights = [ 0 ]
        for i in range(len(self.times)-1, 0, -1):
            deltaT = self.times[i] - self.times[i-1]
            def integrateVal(value, sum, timeSeries):
                sum += deltaT * (value[i-1] + value[i]) / 2
                timeSeries.insert(0, sum)
                return sum

            self.initialOxidizerWeight = integrateVal(oxFlowRate, self.initialOxidizerWeight, self.oxWeights)
            self.initialFuelWeight = integrateVal(fuelFlowRate, self.initialFuelWeight, self.fuelWeights)

    #### Operational Functions ####
    def getInertia(self, time, state):
        timeSinceIgnition = max(0, time - self.ignitionTime)

        oxInertia = self._getOxInertia(timeSinceIgnition)
        fuelInertia = self._getFuelInertia(timeSinceIgnition)
        
        return oxInertia + fuelInertia

    def getAppliedForce(self, state, time, environment, CG):
        #TODO: Model "thrust damping" - where gases moving quickly in the engine act to damp out rotation about the x and y axes
        #TODO: Thrust vs altitude compensation
        timeSinceIgnition = max(0, time - self.ignitionTime)
        
        # Determine thrust magnitude
        if timeSinceIgnition < 0 or timeSinceIgnition > self.times[-1]:
            thrustMagnitude = 0
        else:
            thrustMagnitude = linInterp(self.times, self.thrustLevels, timeSinceIgnition)
        
        # Create Vector
        thrust = Vector(0,0, thrustMagnitude)
        return ForceMomentSystem(thrust)

    def updateIgnitionTime(self, ignitionTime, fakeValue=False):
        self.ignitionTime = ignitionTime
        if not fakeValue:
            self.rocket.engineShutOffTime = max(self.rocket.engineShutOffTime, self.ignitionTime + self.times[-1])
            self.stage.engineShutOffTime = self.ignitionTime + self.times[-1]

    def getTotalImpulse(self):
        # Integrate the thrust - assume linear interpolations between points given -> midpoint rule integrates this perfectly
        totalImpulse = 0
        for i in range(1, len(self.times)):
            deltaT = self.times[i] - self.times[i-1]
            totalImpulse += deltaT * (self.thrustLevels[i-1] + self.thrustLevels[i]) / 2
        
        return totalImpulse

    def _getMass(self, timeSinceIgnition):
        return self.OxWeight(timeSinceIgnition) + self.FuelWeight(timeSinceIgnition)

    def _getOxInertia(self, timeSinceIgnition):
        if self.initialOxidizerWeight == 0:
            return Inertia(Vector(0,0,0), Vector(0,0,0), 0)

        oxWeight = linInterp(self.times, self.oxWeights, timeSinceIgnition)
        
        # Find fraction of oxidizer burned
        oxBurnedFrac = 1 - (oxWeight/self.initialOxidizerWeight)
        
        #Linearly interpolate CG location based on fraction of oxidizer burned
        oxCG_Z = self.initOxCG_Z*(1 - oxBurnedFrac) + self.finalOxCG_Z*oxBurnedFrac
        #TODO: Allow motor(s) to be defined off-axis
        oxCG = Vector(0,0,oxCG_Z)

        # Get MOI
        oxMOI = linInterp(self.times, self.oxMOIs, timeSinceIgnition)
        
        return Inertia(oxMOI, oxCG, oxWeight)

    def _getFuelInertia(self, timeSinceIgnition):
        if self.initialFuelWeight == 0:
            return Inertia(Vector(0,0,0), Vector(0,0,0), 0)

        #See comments in _getOxInertia()
        fuelWeight = linInterp(self.times, self.fuelWeights, timeSinceIgnition)

        fuelBurnedFrac = 1 - (fuelWeight / self.initialFuelWeight)

        fuelCG_Z = self.initFuelCG_Z*(1 - fuelBurnedFrac) + self.finalFuelCG_Z*fuelBurnedFrac
        fuelCG = Vector(0,0,fuelCG_Z)

        fuelMOI = linInterp(self.times, self.fuelMOIs, timeSinceIgnition)

        return Inertia(fuelMOI, fuelCG, fuelWeight)

Ancestors

RocketComponent
abc.ABC

Methods

def getAppliedForce(self, state, time, environment, CG)

Expand source code

def getAppliedForce(self, state, time, environment, CG):
    #TODO: Model "thrust damping" - where gases moving quickly in the engine act to damp out rotation about the x and y axes
    #TODO: Thrust vs altitude compensation
    timeSinceIgnition = max(0, time - self.ignitionTime)
    
    # Determine thrust magnitude
    if timeSinceIgnition < 0 or timeSinceIgnition > self.times[-1]:
        thrustMagnitude = 0
    else:
        thrustMagnitude = linInterp(self.times, self.thrustLevels, timeSinceIgnition)
    
    # Create Vector
    thrust = Vector(0,0, thrustMagnitude)
    return ForceMomentSystem(thrust)

def getInertia(self, time, state)

Expand source code

def getInertia(self, time, state):
    timeSinceIgnition = max(0, time - self.ignitionTime)

    oxInertia = self._getOxInertia(timeSinceIgnition)
    fuelInertia = self._getFuelInertia(timeSinceIgnition)
    
    return oxInertia + fuelInertia

def getTotalImpulse(self)

Expand source code

def getTotalImpulse(self):
    # Integrate the thrust - assume linear interpolations between points given -> midpoint rule integrates this perfectly
    totalImpulse = 0
    for i in range(1, len(self.times)):
        deltaT = self.times[i] - self.times[i-1]
        totalImpulse += deltaT * (self.thrustLevels[i-1] + self.thrustLevels[i]) / 2
    
    return totalImpulse

def updateIgnitionTime(self, ignitionTime, fakeValue=False)

Expand source code

def updateIgnitionTime(self, ignitionTime, fakeValue=False):
    self.ignitionTime = ignitionTime
    if not fakeValue:
        self.rocket.engineShutOffTime = max(self.rocket.engineShutOffTime, self.ignitionTime + self.times[-1])
        self.stage.engineShutOffTime = self.ignitionTime + self.times[-1]

class Transition (*args)

Models a conical diameter transition (growing or shrinking)

Expand source code

class Transition(FixedMass, BodyComponent):
    ''' Models a conical diameter transition (growing or shrinking) '''

    def __init__(self, *args):
        '''
            Two possible sets of inputs:  
            1. Initialization as a regular, dictionary-defined rocket component:  
                * args = (componentDictReader, rocket, stage)  
                * Expected classes: (`MAPLEAF.IO.SubDictReader`, `MAPLEAF.Rocket.Rocket`, `MAPLEAF.Rocket.Stage`)  
            2. Manual initialization:  
                * args = (startDiameter, endDiameter, length, position, inertia, rocket, stage, name, surfaceRoughness)  
                * Expected classes: (float, float, float, `MAPLEAF.Motion.Vector`, `MAPLEAF.Motion.Inertia`, `MAPLEAF.Rocket.Rocket`, `MAPLEAF.Rocket.Stage`, str, float)  
        '''
        if len(args) == 3:
            # Classic initialization from componentDictReader
            componentDictReader, rocket, stage = args
            FixedMass.__init__(self, componentDictReader, rocket, stage)

            self.length = componentDictReader.getFloat("length")
            self.startDiameter = componentDictReader.getFloat("startDiameter")
            self.endDiameter = componentDictReader.getFloat("endDiameter")
            self.surfaceRoughness = componentDictReader.tryGetFloat("surfaceRoughness", defaultValue=self.rocket.surfaceRoughness)
        else:
            # Initialization from parameters passed in
            self.startDiameter, self.endDiameter, self.length, self.position, self.inertia, self.rocket, self.stage, self.name, self.surfaceRoughness = args

        self._precomputeGeometry()

    def _precomputeGeometry(self):
        # Calculate basic areas
        self.topArea = self.startDiameter**2 * math.pi/4
        self.bottomArea = self.endDiameter**2 * math.pi/4
        self.frontalArea = abs(self.topArea - self.bottomArea)

        # Calculate surface area, volume, CP
        if self.length == 0:
            self.wettedArea = 0.0
            self.volume = 0.0
            self.CPLocation = self.position
        else:
            maxDiameter = max(self.startDiameter, self.endDiameter)
            minDiameter = min(self.startDiameter, self.endDiameter)

            hypotenuseSlope = ((maxDiameter - minDiameter) / self.length)
            # Boat tail is a truncated cone
            # Compute surface area of non-truncated cone starting at top of boattail
            baseConeHeight = maxDiameter / hypotenuseSlope
            baseHypotenuse = math.sqrt((maxDiameter/2)**2 + baseConeHeight**2)
            fullConeSurfaceArea = math.pi * (maxDiameter/2) * baseHypotenuse
            # Compute surface area of non-truncated cone starting at bottom of boattail
            tipConeHeight = minDiameter / hypotenuseSlope
            tipHypotenuse = math.sqrt((minDiameter/2)**2 + tipConeHeight**2)
            tipConeSurfaceArea = math.pi * (minDiameter/2) * tipHypotenuse
            # Surface area is the difference b/w the two
            self.wettedArea = fullConeSurfaceArea - tipConeSurfaceArea

            # Volume calculation uses same method as surface area calculation above
            fullConeVolume = self.topArea * baseConeHeight / 3
            tipConeVolume = self.bottomArea * tipConeHeight / 3
            self.volume = fullConeVolume - tipConeVolume

            # Calculate CP
            self.CPLocation = AeroFunctions.Barrowman_GetCPLocation(self.length, self.topArea, self.bottomArea, self.volume, self.position)

        # Calc aspect ratio
        if self.startDiameter == self.endDiameter:
            aspectRatio = 100 # Avoid division by zero
        else:
            aspectRatio = self.length / abs(self.startDiameter - self.endDiameter)

        # Precompute drag properties
        if self.endDiameter <= self.startDiameter:
            # Calculate Cd multiplier based on aspect ratio - Niskanen Eqn 3.88
            if aspectRatio < 1:
                self.CdAdjustmentFactor = 1
            elif aspectRatio < 3:
                self.CdAdjustmentFactor = (3 - aspectRatio) / 2
            else:
                self.CdAdjustmentFactor = 0
        else:
            if self.length == 0:
                coneHalfAngle = math.pi/2
            else:
                coneHalfAngle = math.atan(abs(self.startDiameter - self.endDiameter)/2 / self.length)

            self.coneHalfAngle = coneHalfAngle
            self.SubsonicCdPolyCoeffs = computeSubsonicPolyCoeffs(coneHalfAngle)
            self.TransonicCdPolyCoeffs = computeTransonicPolyCoeffs(coneHalfAngle)

    def plotShape(self):
        import matplotlib.pyplot as plt
        Xvals = []
        Yvals = []

        forePos = self.position.Z
        aftPos = forePos - self.length
        foreRadius = self.startDiameter/2
        aftRadius = self.endDiameter/2

        Xvals.append(forePos)
        Yvals.append(foreRadius) # top right

        Xvals.append(aftPos)
        Yvals.append(aftRadius) # top left

        Xvals.append(aftPos)
        Yvals.append(-aftRadius) # bottom left

        Xvals.append(forePos)
        Yvals.append(-foreRadius) # bottom right

        Xvals.append(forePos)
        Yvals.append(foreRadius) # close in the shape
        plt.plot(Xvals, Yvals, color = 'k')

    def getAppliedForce(self, rocketState, time, environment, CG) -> ForceMomentSystem:
        Mach = AeroParameters.getMachNumber(rocketState, environment)
        Aref = self.rocket.Aref
        
        #### Normal Force ####
        AOA = AeroParameters.getTotalAOA(rocketState, environment)
        CN = AeroFunctions.Barrowman_GetCN(AOA, Aref, self.topArea, self.bottomArea)

        #### Pressure Drag ####
        if self.startDiameter > self.endDiameter:
            # Pressure base drag
            Cd_base = AeroFunctions.getBaseDragCoefficient(Mach)
            Cd_pressure = Cd_base * self.CdAdjustmentFactor

        else:
            # Pressure drag calculated like a nose cone
            if Mach < 1:
                # Niskanen pg. 48 eq. 3.87 - Power law interpolation
                Cd_pressure = self.SubsonicCdPolyCoeffs[0] * Mach**self.SubsonicCdPolyCoeffs[1]

            elif Mach > 1 and Mach < 1.3:
                # Interpolate in transonic region - derived from Niskanen Appendix B, Eqns B.3 - B.6
                Cd_pressure = self.TransonicCdPolyCoeffs[0] + self.TransonicCdPolyCoeffs[1]*Mach +  \
                    self.TransonicCdPolyCoeffs[2]*Mach**2 + self.TransonicCdPolyCoeffs[3]*Mach**3
                
            else:
                Cd_pressure = getSupersonicPressureDragCoeff_Hoerner(self.coneHalfAngle, Mach)
        
        # Make reference are the rocket's, not this objects
        Cd_pressure *= self.frontalArea / Aref 

        #### Skin Friction Drag ####
        if self.wettedArea == 0:
            skinFrictionDragCoefficient = 0
            rollDampingMoment = Vector(0,0,0)
        else:
            skinFrictionDragCoefficient, rollDampingMoment = AeroFunctions.getCylindricalSkinFrictionDragCoefficientAndRollDampingMoment(rocketState, environment, \
                 self.length, Mach, self.surfaceRoughness, self.wettedArea, Aref, self.rocket.finenessRatio, self.rocket.fullyTurbulentBL)

        #### Total Drag ####
        Cd = Cd_pressure + skinFrictionDragCoefficient

        #### Assemble & return final force object ####
        return AeroFunctions.forceFromCdCN(rocketState, environment, Cd, CN, self.CPLocation, Aref, moment=rollDampingMoment)

    def getMaxDiameter(self):
        return max(self.startDiameter, self.endDiameter)

    def getRadius(self, distanceFromTip):
        return (distanceFromTip/self.length * (self.endDiameter - self.startDiameter) + self.startDiameter) / 2

Ancestors

Subclasses

BoatTail

Methods

def getAppliedForce(self, rocketState, time, environment, CG) ‑> ForceMomentSystem

Expand source code

def getAppliedForce(self, rocketState, time, environment, CG) -> ForceMomentSystem:
    Mach = AeroParameters.getMachNumber(rocketState, environment)
    Aref = self.rocket.Aref
    
    #### Normal Force ####
    AOA = AeroParameters.getTotalAOA(rocketState, environment)
    CN = AeroFunctions.Barrowman_GetCN(AOA, Aref, self.topArea, self.bottomArea)

    #### Pressure Drag ####
    if self.startDiameter > self.endDiameter:
        # Pressure base drag
        Cd_base = AeroFunctions.getBaseDragCoefficient(Mach)
        Cd_pressure = Cd_base * self.CdAdjustmentFactor

    else:
        # Pressure drag calculated like a nose cone
        if Mach < 1:
            # Niskanen pg. 48 eq. 3.87 - Power law interpolation
            Cd_pressure = self.SubsonicCdPolyCoeffs[0] * Mach**self.SubsonicCdPolyCoeffs[1]

        elif Mach > 1 and Mach < 1.3:
            # Interpolate in transonic region - derived from Niskanen Appendix B, Eqns B.3 - B.6
            Cd_pressure = self.TransonicCdPolyCoeffs[0] + self.TransonicCdPolyCoeffs[1]*Mach +  \
                self.TransonicCdPolyCoeffs[2]*Mach**2 + self.TransonicCdPolyCoeffs[3]*Mach**3
            
        else:
            Cd_pressure = getSupersonicPressureDragCoeff_Hoerner(self.coneHalfAngle, Mach)
    
    # Make reference are the rocket's, not this objects
    Cd_pressure *= self.frontalArea / Aref 

    #### Skin Friction Drag ####
    if self.wettedArea == 0:
        skinFrictionDragCoefficient = 0
        rollDampingMoment = Vector(0,0,0)
    else:
        skinFrictionDragCoefficient, rollDampingMoment = AeroFunctions.getCylindricalSkinFrictionDragCoefficientAndRollDampingMoment(rocketState, environment, \
             self.length, Mach, self.surfaceRoughness, self.wettedArea, Aref, self.rocket.finenessRatio, self.rocket.fullyTurbulentBL)

    #### Total Drag ####
    Cd = Cd_pressure + skinFrictionDragCoefficient

    #### Assemble & return final force object ####
    return AeroFunctions.forceFromCdCN(rocketState, environment, Cd, CN, self.CPLocation, Aref, moment=rollDampingMoment)

def plotShape(self)

Expand source code

def plotShape(self):
    import matplotlib.pyplot as plt
    Xvals = []
    Yvals = []

    forePos = self.position.Z
    aftPos = forePos - self.length
    foreRadius = self.startDiameter/2
    aftRadius = self.endDiameter/2

    Xvals.append(forePos)
    Yvals.append(foreRadius) # top right

    Xvals.append(aftPos)
    Yvals.append(aftRadius) # top left

    Xvals.append(aftPos)
    Yvals.append(-aftRadius) # bottom left

    Xvals.append(forePos)
    Yvals.append(-foreRadius) # bottom right

    Xvals.append(forePos)
    Yvals.append(foreRadius) # close in the shape
    plt.plot(Xvals, Yvals, color = 'k')

Inherited members

BodyComponent:
- getBottomInterfaceLocation
- getMaxDiameter
- getRadius
- getTopInterfaceLocation