# ISC License
#
# Copyright (c) 2016-2017, 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 os
import numpy as np
import pytest
from Basilisk.architecture import messaging
from Basilisk.simulation import msisAtmosphere
# 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 # general support file with common unit test functions
filename = inspect.getframeinfo(inspect.currentframe()).filename
path = os.path.dirname(os.path.abspath(filename))
bskName = 'Basilisk'
splitPath = path.split(bskName)
# 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.
[docs]@pytest.mark.parametrize("orbitType", ["LPO", "LTO"])
@pytest.mark.parametrize("setEpoch", ["Default", "Direct", "Msg"])
# provide a unique test method name, starting with test_
def test_scenarioMsisAtmosphereOrbit(show_plots, orbitType, setEpoch):
"""This function is called by the py.test environment."""
# each test method requires a single assert method to be called
showVal = False
[testResults, testMessage] = run(showVal, orbitType, setEpoch)
assert testResults < 1, testMessage
[docs]def run(show_plots, orbitCase, setEpoch):
"""Call this routine directly to run the script."""
testFailCount = 0 # zero unit test result counter
testMessages = [] # create empty array to store test log messages
#
# From here on there scenario python code is found. Above this line the code is to setup a
# unitTest environment. The above code is not critical if learning how to code BSK.
#
# 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(10.)
dynProcess.addTask(scSim.CreateNewTask(simTaskName, simulationTimeStep))
# Initialize new atmosphere and drag model, add them to task
newAtmo = msisAtmosphere.MsisAtmosphere()
atmoTaskName = "atmosphere"
newAtmo.ModelTag = "MsisAtmo"
if setEpoch == "Msg":
epochMsg = unitTestSupport.timeStringToGregorianUTCMsg('2019 Jan 01 00:00:00.00 (UTC)')
newAtmo.epochInMsg.subscribeTo(epochMsg)
# setting epoch day of year info deliberately to a false value. The epoch msg info should be used
newAtmo.epochDoy = 10
elif setEpoch == "Direct":
newAtmo.epochDoy = 1 # setting epoch day of year info directly
dynProcess.addTask(scSim.CreateNewTask(atmoTaskName, simulationTimeStep))
scSim.AddModelToTask(atmoTaskName, newAtmo)
# initialize spacecraft object and set properties
scObject = spacecraft.Spacecraft()
scObject.ModelTag = "spacecraftBody"
# add spacecraft object to the simulation process
scSim.AddModelToTask(simTaskName, scObject)
# clear prior gravitational body and SPICE setup definitions
gravFactory = simIncludeGravBody.gravBodyFactory()
newAtmo.addSpacecraftToModel(scObject.scStateOutMsg)
planet = gravFactory.createEarth()
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()
r_eq = planet.radEquator
if orbitCase == "LPO":
orbAltMin = 100.0*1000.0
orbAltMax = orbAltMin
elif orbitCase == "LTO":
orbAltMin = 100.*1000.0
orbAltMax = 100.0 * 1000.0
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(0.002*P)
#
# Setup data logging before the simulation is initialized
#
sw_msg_names = [
"ap_24_0", "ap_3_0", "ap_3_-3", "ap_3_-6", "ap_3_-9",
"ap_3_-12", "ap_3_-15", "ap_3_-18", "ap_3_-21", "ap_3_-24",
"ap_3_-27", "ap_3_-30", "ap_3_-33", "ap_3_-36", "ap_3_-39",
"ap_3_-42", "ap_3_-45", "ap_3_-48", "ap_3_-51", "ap_3_-54",
"ap_3_-57", "f107_1944_0", "f107_24_-24"
]
swMsgList = []
for c in range(len(sw_msg_names)):
swMsgData = messaging.SwDataMsgPayload()
swMsgData.dataValue = 0
swMsgList.append(messaging.SwDataMsg().write(swMsgData))
newAtmo.swDataInMsgs[c].subscribeTo(swMsgList[-1])
dataLog = scObject.scStateOutMsg.recorder()
denLog = newAtmo.envOutMsgs[0].recorder()
scSim.AddModelToTask(simTaskName, dataLog)
scSim.AddModelToTask(simTaskName, denLog)
#
# initialize Spacecraft States with the 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
densData = denLog.neutralDensity
tempData = denLog.localTemp
np.set_printoptions(precision=16)
# Compare to expected values
refAtmoData = np.loadtxt(path + '/truthOutputs.txt')
accuracy = 1e-8
unitTestSupport.writeTeXSnippet("unitTestToleranceValue", str(accuracy), path)
# Test atmospheric density calculation; note that refAtmoData is in g/cm^3,
# and must be adjusted by a factor of 1e-3 to match kg/m^3
print(densData[-1])
print(refAtmoData[5]*1000)
if np.testing.assert_allclose(densData[-1], refAtmoData[5]*1000., atol=accuracy):
testFailCount += 1
testMessages.append("FAILED: NRLMSISE-00 failed density unit test with a value difference of "+str(densData[0]-refAtmoData[5]*1000))
print(tempData[-1])
print(refAtmoData[-1])
if np.testing.assert_allclose(tempData[-1], refAtmoData[-1], atol=accuracy):
testFailCount += 1
testMessages.append(
"FAILED: NRLMSISE-00 failed temperature unit test with a value difference of "+str(tempData[-1]-refAtmoData[-1]))
snippentName = "unitTestPassFail" + str(orbitCase) + str(setEpoch)
if testFailCount == 0:
colorText = 'ForestGreen'
print("PASSED: " + newAtmo.ModelTag)
passedText = '\\textcolor{' + colorText + '}{' + "PASSED" + '}'
else:
colorText = 'Red'
print("Failed: " + newAtmo.ModelTag)
passedText = '\\textcolor{' + colorText + '}{' + "Failed" + '}'
print(testMessages)
unitTestSupport.writeTeXSnippet(snippentName, passedText, path)
return [testFailCount, ''.join(testMessages)]
if __name__ == '__main__':
run(True,
"LPO", # orbitCase
"Msg") # setEpoch