#
# ISC License
#
# 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.
#
r"""
Overview
--------
This script sets up a 3-DOF spacecraft which is orbiting Earth. The purpose
is to illustrate how to start and stop the Basilisk simulation to apply
some :math:`\Delta v`'s for simple orbit maneuvers. Read :ref:`scenarioBasicOrbit`
to learn how to setup an orbit simulation.
The script is found in the folder ``basilisk/examples`` and executed by using::
python3 scenarioOrbitManeuver.py
The simulation layout is shown in the following illustration. A single simulation process is created
which contains the spacecraft object. The BSK simulation is run for a fixed period. After stopping, the
states are changed and the simulation is resumed.
.. image:: /_images/static/test_scenarioOrbitManeuver.svg
:align: center
When the simulation completes 2 plots are shown for each case. One plot always shows
the inertial position vector components, while the second plot either shows a plot
of the radius time history (Hohmann maneuver), or the
inclination angle time history plot (Inclination change maneuver).
Illustration of Simulation Results
----------------------------------
The following images illustrate the expected simulation run returns for a range of script configurations.
::
show_plots = True, maneuverCase = 0
In this case a classical Hohmann transfer is being
simulated to go from LEO to reach and stay at GEO. The math behind such maneuvers can be found
in textbooks such as `Analytical Mechanics of Space Systems <http://arc.aiaa.org/doi/book/10.2514/4.102400>`__.
.. image:: /_images/Scenarios/scenarioOrbitManeuver10.svg
:align: center
.. image:: /_images/Scenarios/scenarioOrbitManeuver20.svg
:align: center
::
show_plots = True, maneuverCase = 1
In this case a classical plane change is being
simulated to go rotate the orbit plane first 8 degrees, then another 4 degrees after
orbiting 90 degrees. The math behind such maneuvers can be found
in textbooks such as `Analytical Mechanics of Space Systems
<http://arc.aiaa.org/doi/book/10.2514/4.102400>`__.
.. image:: /_images/Scenarios/scenarioOrbitManeuver11.svg
:align: center
.. image:: /_images/Scenarios/scenarioOrbitManeuver21.svg
:align: center
"""
#
# Basilisk Scenario Script and Integrated Test
#
# Purpose: Integrated test of the spacecraft() and gravity modules illustrating
# how impulsive Delta_v maneuver can be simulated with stopping and starting the
# simulation.
# Author: Hanspeter Schaub
# Creation Date: Nov. 26, 2016
#
import os
import numpy as np
import math
from Basilisk.utilities import SimulationBaseClass
from Basilisk.utilities import unitTestSupport # general support file with common unit test functions
import matplotlib.pyplot as plt
from Basilisk.utilities import macros
from Basilisk.utilities import orbitalMotion
from Basilisk.simulation import spacecraft
from Basilisk.utilities import simIncludeGravBody
from Basilisk.architecture import messaging
# attempt to import vizard
from Basilisk.utilities import vizSupport
# The path to the location of Basilisk
# Used to get the location of supporting data.
from Basilisk import __path__
bskPath = __path__[0]
fileName = os.path.basename(os.path.splitext(__file__)[0])
[docs]def run(show_plots, maneuverCase):
"""
The scenarios can be run with the followings setups parameters:
Args:
show_plots (bool): Determines if the script should display plots
maneuverCase (int):
====== ============================
Int Definition
====== ============================
0 Hohmann maneuver
1 Inclination change maneuver
====== ============================
"""
# 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))
#
# setup the simulation tasks/objects
#
# initialize spacecraft object and set properties
scObject = spacecraft.Spacecraft()
scObject.ModelTag = "spacecraftBody"
# add spacecraft object to the simulation process
scSim.AddModelToTask(simTaskName, scObject)
# setup Gravity Body
gravFactory = simIncludeGravBody.gravBodyFactory()
earth = gravFactory.createEarth()
earth.isCentralBody = True # ensure this is the central gravitational body
# attach gravity model to spacecraft
scObject.gravField.gravBodies = spacecraft.GravBodyVector(list(gravFactory.gravBodies.values()))
#
# setup orbit and simulation time
#
# setup the orbit using classical orbit elements
oe = orbitalMotion.ClassicElements()
rLEO = 7000. * 1000 # meters
rGEO = math.pow(earth.mu / math.pow((2. * np.pi) / (24. * 3600.), 2), 1. / 3.)
oe.a = rLEO
oe.e = 0.0001
oe.i = 0.0 * macros.D2R
oe.Omega = 48.2 * macros.D2R
oe.omega = 347.8 * macros.D2R
oe.f = 85.3 * macros.D2R
rN, vN = orbitalMotion.elem2rv(earth.mu, oe)
scObject.hub.r_CN_NInit = rN # m - r_CN_N
scObject.hub.v_CN_NInit = vN # m - v_CN_N
# set the simulation time
n = np.sqrt(earth.mu / oe.a / oe.a / oe.a)
P = 2. * np.pi / n
simulationTime = macros.sec2nano(0.25 * P)
#
# Setup data logging before the simulation is initialized
#
numDataPoints = 100
samplingTime = unitTestSupport.samplingTime(simulationTime, simulationTimeStep, numDataPoints)
dataRec = scObject.scStateOutMsg.recorder(samplingTime)
scSim.AddModelToTask(simTaskName, dataRec)
# if this scenario is to interface with the BSK Viz, uncomment the following lines
vizSupport.enableUnityVisualization(scSim, simTaskName, scObject
# , saveFile=fileName
)
#
# initialize Simulation
#
scSim.InitializeSimulation()
#
# get access to dynManager translational states for future access to the states
#
posRef = scObject.dynManager.getStateObject("hubPosition")
velRef = scObject.dynManager.getStateObject("hubVelocity")
# The dynamics simulation is setup using a Spacecraft() module with the Earth's
# gravity module attached. Note that the rotational motion simulation is turned off to leave
# pure 3-DOF translation motion simulation. After running the simulation for 1/4 of a period
# the simulation is stopped to apply impulsive changes to the inertial velocity vector.
scSim.ConfigureStopTime(simulationTime)
scSim.ExecuteSimulation()
# Next, the state manager objects are called to retrieve the latest inertial position and
# velocity vector components:
rVt = unitTestSupport.EigenVector3d2np(posRef.getState())
vVt = unitTestSupport.EigenVector3d2np(velRef.getState())
# compute maneuver Delta_v's
if maneuverCase == 1:
# inclination change
Delta_i = 8.0 * macros.D2R
rHat = rVt / np.linalg.norm(rVt)
hHat = np.cross(rVt, vVt)
hHat = hHat / np.linalg.norm(hHat)
vHat = np.cross(hHat, rHat)
v0 = np.dot(vHat, vVt)
vVt = vVt - (1. - np.cos(Delta_i)) * v0 * vHat + np.sin(Delta_i) * v0 * hHat
# After computing the maneuver specific Delta_v's, the state managers velocity is updated through
velRef.setState(unitTestSupport.np2EigenVectorXd(vVt))
T2 = macros.sec2nano(P * 0.25)
else:
# Hohmann Transfer to GEO
v0 = np.linalg.norm(vVt)
r0 = np.linalg.norm(rVt)
at = (r0 + rGEO) * .5
v0p = np.sqrt(earth.mu / at * rGEO / r0)
n1 = np.sqrt(earth.mu / at / at / at)
T2 = macros.sec2nano((np.pi) / n1)
vHat = vVt / v0
vVt = vVt + vHat * (v0p - v0)
# After computing the maneuver specific Delta_v's, the state managers velocity is updated through
velRef.setState(unitTestSupport.np2EigenVectorXd(vVt))
# To start up the simulation again, note that the total simulation time must be provided,
# not just the next incremental simulation time.
scSim.ConfigureStopTime(simulationTime + T2)
scSim.ExecuteSimulation()
# This process is then repeated for the second maneuver.
# get the current spacecraft states
rVt = unitTestSupport.EigenVector3d2np(posRef.getState())
vVt = unitTestSupport.EigenVector3d2np(velRef.getState())
# compute maneuver Delta_v's
if maneuverCase == 1:
# inclination change
Delta_i = 4.0 * macros.D2R
rHat = rVt / np.linalg.norm(rVt)
hHat = np.cross(rVt, vVt)
hHat = hHat / np.linalg.norm(hHat)
vHat = np.cross(hHat, rHat)
v0 = np.dot(vHat, vVt)
vVt = vVt - (1. - np.cos(Delta_i)) * v0 * vHat + np.sin(Delta_i) * v0 * hHat
velRef.setState(unitTestSupport.np2EigenVectorXd(vVt))
T3 = macros.sec2nano(P * 0.25)
else:
# Hohmann Transfer to GEO
v1 = np.linalg.norm(vVt)
v1p = np.sqrt(earth.mu / rGEO)
n1 = np.sqrt(earth.mu / rGEO / rGEO / rGEO)
T3 = macros.sec2nano(0.25 * (np.pi) / n1)
vHat = vVt / v1
vVt = vVt + vHat * (v1p - v1)
velRef.setState(unitTestSupport.np2EigenVectorXd(vVt))
# run simulation for 3rd chunk
scSim.ConfigureStopTime(simulationTime + T2 + T3)
scSim.ExecuteSimulation()
#
# retrieve the logged data
#
posData = dataRec.r_BN_N
velData = dataRec.v_BN_N
np.set_printoptions(precision=16)
#
# plot the results
#
# draw the inertial position vector components
plt.close("all") # clears out plots from earlier test runs
plt.figure(1)
fig = plt.gcf()
ax = fig.gca()
ax.ticklabel_format(useOffset=False, style='plain')
for idx in range(3):
plt.plot(dataRec.times() * macros.NANO2HOUR, posData[:, idx] / 1000.,
color=unitTestSupport.getLineColor(idx, 3),
label='$r_{BN,' + str(idx) + '}$')
plt.legend(loc='lower right')
plt.xlabel('Time [h]')
plt.ylabel('Inertial Position [km]')
figureList = {}
pltName = fileName + "1" + str(int(maneuverCase))
figureList[pltName] = plt.figure(1)
if maneuverCase == 1:
# show inclination angle
plt.figure(2)
fig = plt.gcf()
ax = fig.gca()
ax.ticklabel_format(useOffset=False, style='plain')
iData = []
for idx in range(0, len(posData)):
oeData = orbitalMotion.rv2elem(earth.mu, posData[idx], velData[idx])
iData.append(oeData.i * macros.R2D)
plt.plot(dataRec.times() * macros.NANO2HOUR, np.ones(len(posData[:, 0])) * 8.93845, '--', color='#444444'
)
plt.plot(dataRec.times() * macros.NANO2HOUR, iData, color='#aa0000'
)
plt.ylim([-1, 10])
plt.xlabel('Time [h]')
plt.ylabel('Inclination [deg]')
else:
# show SMA
plt.figure(2)
fig = plt.gcf()
ax = fig.gca()
ax.ticklabel_format(useOffset=False, style='plain')
rData = []
for idx in range(0, len(posData)):
oeData = orbitalMotion.rv2elem_parab(earth.mu, posData[idx], velData[idx])
rData.append(oeData.rmag / 1000.)
plt.plot(dataRec.times() * macros.NANO2HOUR, rData, color='#aa0000',
)
plt.xlabel('Time [h]')
plt.ylabel('Radius [km]')
pltName = fileName + "2" + str(int(maneuverCase))
figureList[pltName] = plt.figure(2)
if show_plots:
plt.show()
# close the plots being saved off to avoid over-writing old and new figures
plt.close("all")
# each test method requires a single assert method to be called
# this check below just makes sure no sub-test failures were found
dataPos = posRef.getState()
dataPos = [[0.0, dataPos[0][0], dataPos[1][0], dataPos[2][0]]]
return figureList
#
# This statement below ensures that the unit test scrip can be run as a
# stand-along python script
#
if __name__ == "__main__":
run(
True, # show_plots
0 # Maneuver Case (0 - Hohmann, 1 - Inclination)
)