Source code for test_atmoDrag


# Copyright (c) 2016, Autonomous Vehicle Systems Lab, University of Colorado at Boulder
#
# Permission to use, copy, modify, and/or distribute this software for any
# purpose with or without fee is hereby granted, provided that the above
# copyright notice and this permission notice appear in all copies.
#
# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
# WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
# MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
# ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
# WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
# ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
# OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.


#
# Basilisk Scenario Script and Integrated Test
#
# Purpose:  Test the validity of a simple exponential atmosphere model.
# Author:   Andrew Harris
# Creation Date:  Jan 18, 2017
#

import inspect
import math
import os

import matplotlib.pyplot as plt
import numpy as np
# print dir(exponentialAtmosphere)
from Basilisk.simulation import dragDynamicEffector
from Basilisk.simulation import exponentialAtmosphere, simpleNav
# import simulation related support
from Basilisk.simulation import spacecraft
# import general simulation support files
from Basilisk.utilities import SimulationBaseClass
from Basilisk.utilities import macros
from Basilisk.utilities import orbitalMotion
from Basilisk.utilities import simIncludeGravBody
from Basilisk.utilities import unitTestSupport, RigidBodyKinematics

filename = inspect.getframeinfo(inspect.currentframe()).filename
path = os.path.dirname(os.path.abspath(filename))


# uncomment this line is this test is to be skipped in the global unit test run, adjust message as needed
# @pytest.mark.skipif(conditionstring)
# uncomment this line if this test has an expected failure, adjust message as needed
# @pytest.mark.xfail(True, reason="Previously set sim parameters are not consistent with new formulation\n")

# The following 'parametrize' function decorator provides the parameters and expected results for each
#   of the multiple test runs for this test.

# provide a unique test method name, starting with test_
[docs] def test_scenarioDragOrbit(): """This function is called by the py.test environment.""" # each test method requires a single assert method to be called earthCase = "Earth" marsCase = "Mars" orb1 = "LPO" orb2 = "LTO" showVal = False testResults = [] testMessage = [] [leoResults, leoMessages] = run( showVal, orb1, earthCase) #[gtoResults, gtoMessages] = run( # showVal, orb2, earthCase) #[lmoResults, lmoMessages] = run( #showVal, orb1, marsCase) #[mtoResults, mtoMessages] = run( #showVal, orb2, marsCase) testResults = leoResults#+gtoResults#+lmoResults+mtoResults testMessage.append(leoMessages) #testMessage.append(gtoMessages) ##testMessage.append(lmoMessages) #testMessage.append(mtoMessages) assert testResults < 1, testMessage
def expAtmoComp(alt, baseDens, scaleHeight): dens = baseDens * math.exp(-alt/scaleHeight) return dens def cannonballDragComp(dragCoeff, dens, area, vel, att): dragDir_N = -vel / np.linalg.norm(vel) dcm_BN = RigidBodyKinematics.MRP2C(att) dragDir_B = dcm_BN.dot(dragDir_N) dragForce = 0.5 * dragCoeff * dens * area * np.linalg.norm(vel)**2.0 * dragDir_B return dragForce
[docs] def run(show_plots, orbitCase, planetCase): """Call this routine directly to run the tutorial scenario.""" testFailCount = 0 # zero unit test result counter testMessages = [] # create empty array to store test log messages # Create simulation variable names simTaskName = "simTask" simProcessName = "simProcess" # Create a sim module as an empty container scSim = SimulationBaseClass.SimBaseClass() # create the simulation process dynProcess = scSim.CreateNewProcess(simProcessName) # create the dynamics task and specify the integration update time simulationTimeStep = macros.sec2nano(1.0) dynProcess.addTask(scSim.CreateNewTask(simTaskName, simulationTimeStep)) # Initialize new atmosphere and drag model, add them to task newAtmo = exponentialAtmosphere.ExponentialAtmosphere() atmoTaskName = "atmosphere" newAtmo.ModelTag = "ExpAtmo" projArea = 10.0 # Set drag area in m^2 dragCoeff = 2.2 # Set drag ceofficient dragEffector = dragDynamicEffector.DragDynamicEffector() dragEffector.ModelTag = "DragEff" dragEffectorTaskName = "drag" dragEffector.coreParams.projectedArea = projArea dragEffector.coreParams.dragCoeff = dragCoeff dragEffector.coreParams.comOffset = [1., 0., 0.] dynProcess.addTask(scSim.CreateNewTask(atmoTaskName, simulationTimeStep)) dynProcess.addTask(scSim.CreateNewTask(dragEffectorTaskName, simulationTimeStep)) scSim.AddModelToTask(atmoTaskName, newAtmo) # # setup the simulation tasks/objects # # initialize spacecraft object and set properties scObject = spacecraft.Spacecraft() scObject.ModelTag = "spacecraftBody" newAtmo.addSpacecraftToModel(scObject.scStateOutMsg) simpleNavObj = simpleNav.SimpleNav() scSim.AddModelToTask(simTaskName, simpleNavObj) simpleNavObj.scStateInMsg.subscribeTo(scObject.scStateOutMsg) scObject.addDynamicEffector(dragEffector) # add spacecraft object to the simulation process scSim.AddModelToTask(simTaskName, scObject) scSim.AddModelToTask(dragEffectorTaskName, dragEffector) # clear prior gravitational body and SPICE setup definitions gravFactory = simIncludeGravBody.gravBodyFactory() dragEffector.atmoDensInMsg.subscribeTo(newAtmo.envOutMsgs[0]) if planetCase == "Earth": planet = gravFactory.createEarth() elif planetCase == "Mars": planet = gravFactory.createMars() planet.isCentralBody = True # ensure this is the central gravitational body mu = planet.mu # attach gravity model to spacecraft scObject.gravField.gravBodies = spacecraft.GravBodyVector(list(gravFactory.gravBodies.values())) # # setup orbit and simulation time oe = orbitalMotion.ClassicElements() if planetCase == "Earth": r_eq = 6371*1000.0 refBaseDens = 1.217 refScaleHeight = 8500.0 elif planetCase == "Mars": refBaseDens = 0.020 refScaleHeight = 11100.0 r_eq = 3389.5 * 1000.0 else: return 1, "Test failed- did not initialize planets." if orbitCase == "LPO": orbAltMin = 300.0*1000.0 orbAltMax = orbAltMin elif orbitCase == "LTO": orbAltMin = 300*1000.0 orbAltMax = 800.0 * 1000.0 newAtmo.planetRadius = r_eq newAtmo.scaleHeight = refScaleHeight newAtmo.baseDensity = refBaseDens rMin = r_eq + orbAltMin rMax = r_eq + orbAltMax oe.a = (rMin+rMax)/2.0 oe.e = 1.0 - rMin/oe.a oe.i = 0.0*macros.D2R oe.Omega = 0.0*macros.D2R oe.omega = 0.0*macros.D2R oe.f = 0.0*macros.D2R rN, vN = orbitalMotion.elem2rv(mu, oe) oe = orbitalMotion.rv2elem(mu, rN, vN) # this stores consistent initial orbit elements # with circular or equatorial orbit, some angles are # arbitrary # set the simulation time n = np.sqrt(mu/oe.a/oe.a/oe.a) P = 2.*np.pi/n simulationTime = macros.sec2nano(1*P) # # Setup data logging before the simulation is initialized # numDataPoints = 100 samplingTime = unitTestSupport.samplingTime(simulationTime, simulationTimeStep, numDataPoints) dataLog = scObject.scStateOutMsg.recorder(samplingTime) dataNewAtmoLog = newAtmo.envOutMsgs[0].recorder(samplingTime) scSim.AddModelToTask(simTaskName, dataLog) scSim.AddModelToTask(simTaskName, dataNewAtmoLog) dragEffectorLog = dragEffector.logger("forceExternal_B", samplingTime) scSim.AddModelToTask(simTaskName, dragEffectorLog) # # initialize Spacecraft States with initialization variables # scObject.hub.r_CN_NInit = rN # m - r_CN_N scObject.hub.v_CN_NInit = vN # m - v_CN_N # # initialize Simulation # scSim.InitializeSimulation() # # configure a simulation stop time and execute the simulation run # scSim.ConfigureStopTime(simulationTime) scSim.ExecuteSimulation() # # retrieve the logged data # posData = dataLog.r_BN_N velData = dataLog.v_BN_N attData = dataLog.sigma_BN dragForce = dragEffectorLog.forceExternal_B densData = dataNewAtmoLog.neutralDensity np.set_printoptions(precision=16) # Compare to expected values endInd = dragForce.shape[0] refDragForce = np.zeros([endInd,3]) refDensData = np.zeros([endInd,1]) accuracy = 1e-13 # print planetCase # print orbitCase for ind in range(0, endInd-1): # print "Position data:", posData[ind,1:] # print "Velocity data:", velData[ind,1:] # print "Density data:", densData[ind,1] refDragForce[ind] = cannonballDragComp(dragCoeff,densData[ind],projArea,velData[ind], attData[ind]) # print "Reference drag data:", refDragForce[ind,:] # print "Drag Data:", dragForce[ind,:] # print "" # check a vector values for ind in range(1,endInd-1): if not unitTestSupport.isArrayEqual(dragForce[ind,:], refDragForce[ind],3,accuracy): testFailCount += 1 testMessages.append( "FAILED: DragEffector failed force unit test with a value difference of " + str(np.linalg.norm(dragForce[ind,:]-refDragForce[ind]))) # # plot the results # if show_plots: plt.close("all") # clears out plots from earlier test runs # draw the inertial position vector components plt.figure(1) fig = plt.gcf() ax = fig.gca() ax.ticklabel_format(useOffset=False, style='plain') for idx in range(0,3): plt.plot(dataLog.times()*macros.NANO2SEC/P, posData[:, idx]/1000., color=unitTestSupport.getLineColor(idx,3), label='$r_{BN,'+str(idx)+'}$') plt.legend(loc='lower right') plt.xlabel('Time [orbits]') plt.ylabel('Inertial Position [km]') # draw orbit in perifocal frame b = oe.a*np.sqrt(1-oe.e*oe.e) p = oe.a*(1-oe.e*oe.e) plt.figure(2,figsize=np.array((1.0, b/oe.a))*4.75,dpi=100) plt.axis(np.array([-oe.rApoap, oe.rPeriap, -b, b])/1000*1.25) # draw the planet fig = plt.gcf() ax = fig.gca() planetColor= '#008800' planetRadius = planet.radEquator/1000 ax.add_artist(plt.Circle((0, 0), planetRadius, color=planetColor)) # draw the actual orbit rData=[] fData=[] for idx in range(0,len(posData)): oeData = orbitalMotion.rv2elem(mu,posData[idx,0:3],velData[idx,0:3]) rData.append(oeData.rmag) fData.append(oeData.f + oeData.omega - oe.omega) plt.plot(rData*np.cos(fData)/1000, rData*np.sin(fData)/1000 ,color='#aa0000' ,linewidth = 3.0 ) # draw the full osculating orbit from the initial conditions fData = np.linspace(0,2*np.pi,100) rData = [] for idx in range(0,len(fData)): rData.append(p/(1+oe.e*np.cos(fData[idx]))) plt.plot(rData*np.cos(fData)/1000, rData*np.sin(fData)/1000 ,'--' , color='#555555' ) plt.xlabel('$i_e$ Cord. [km]') plt.ylabel('$i_p$ Cord. [km]') plt.grid() plt.figure() fig = plt.gcf() ax = fig.gca() ax.ticklabel_format(useOffset=False, style='plain') smaData = [] for idx in range(0, len(posData)): oeData = orbitalMotion.rv2elem(mu, posData[idx, 0:3], velData[idx, 0:3]) smaData.append(oeData.a/1000.) plt.plot(posData[:, 0]*macros.NANO2SEC/P, smaData ,color='#aa0000', ) plt.xlabel('Time [orbits]') plt.ylabel('SMA [km]') plt.figure() fig = plt.gcf() ax = fig.gca() ax.ticklabel_format(useOffset=False, style='sci') plt.plot( dataNewAtmoLog.times()*macros.NANO2SEC, densData) plt.title('Density Data vs. Time') plt.xlabel('Time') plt.ylabel('Density in kg/m^3') plt.show() plt.close("all") if testFailCount == 0: print("PASSED: " + dragEffector.ModelTag) else: print("Failed: " + dragEffector.ModelTag) return testFailCount, testMessages
# close the plots being saved off to avoid over-writing old and new figures if __name__ == '__main__': run(True,"LPO","Earth")