#
# ISC License
#
# Copyright (c) 2022, 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
--------
Demonstrates the interaction between two charged spacecraft in a leader/follower configuration, the effect electrostatic forces/torques have on the separation of the two spacecrafts and how to visualize the simulation data in :ref:`Vizard <vizard>`. The scenario demonstrates how to use :ref:`msmForceTorque` to calculate the electrostatic forces using the Multi-Sphere Method. Each spacecraft is represented by multiple spheres each with a designated location and radius. The locations and radii data is stored in ``GOESR_bus_80_sphs.csv``, however any appropriate csv file can be used. Both spacecraft have a negative charged potential. The purpose of this script is to show how to set up the Multi-Sphere Method for charged spacecraft and apply the external forces/torques to the spacecrafts as well as to show how to store the Basilisk simulation data to be able to visualize both satellite's motions within the :ref:`Vizard <vizard>` application.
The script is found in the folder ``basilisk/examples`` and executed by using::
python3 scenarioTwoChargedSC.py
The simulation layout is shown in the following illustration. A single simulation process is created
which contains both the leader spacecraft and the follower object as well as the multisphere model for each of the spacecrafts.
.. image:: /_images/static/test_scenarioTwoChargedSC.svg
:align: center
When the simulation completes, several plots are shown for the separation distance of the two satellites, the relative orbit of the follower spacecraft around the leader spacecraft, and the Multisphere Model representation for both spacecraft with a color bar that denotes charge of the spheres.
Illustration of Simulation Results
----------------------------------
::
show_plots = True
.. image:: /_images/Scenarios/scenarioTwoChargedSC1.svg
:align: center
.. image:: /_images/Scenarios/scenarioTwoChargedSC2.svg
:align: center
.. image:: /_images/Scenarios/scenarioTwoChargedSC3.svg
:align: center
.. image:: /_images/Scenarios/scenarioTwoChargedSC4.svg
:align: center
"""
#
# Basilisk Scenario Script and Integrated Test
#
# Purpose: Basic simulation showing two charged spacecraft interacting based on the Multisphere Sphere Method.
# Author: James Walker
# Creation Date: January 19, 2022
#
import copy
import csv
import math
import os
import matplotlib.pyplot as plt
import numpy as np
from Basilisk.architecture import messaging
from Basilisk.simulation import spacecraft, extForceTorque, msmForceTorque
from Basilisk.utilities import (SimulationBaseClass, macros,
orbitalMotion, simIncludeGravBody,
unitTestSupport, RigidBodyKinematics, vizSupport, SpherePlot)
# 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):
"""
The scenarios can be run with the following setup parameters:
Args:
show_plots (bool): Determines if the script should display plots
"""
# Create simulation variable names
dynTaskName = "dynTask"
dynProcessName = "dynProcess"
# Create a sim module as an empty container
scSim = SimulationBaseClass.SimBaseClass()
#
# create the simulation process
#
dynProcess = scSim.CreateNewProcess(dynProcessName, 1)
# create the dynamics task and specify the integration update time
simulationTimeStep = macros.sec2nano(25.0)
dynProcess.addTask(scSim.CreateNewTask(dynTaskName, simulationTimeStep))
#
# setup the simulation tasks/objects
#
# initialize leader spacecraft object and set properties
scObjectLeader = spacecraft.Spacecraft()
scObjectLeader.ModelTag = "Leader"
# initialize follower spacecraft object and set properties
scObjectFollower = spacecraft.Spacecraft()
scObjectFollower.ModelTag = "Follower"
# add spacecraftPlus object to the simulation process
scSim.AddModelToTask(dynTaskName, scObjectLeader)
scSim.AddModelToTask(dynTaskName, scObjectFollower)
# clear prior gravitational body and SPICE setup definitions
gravFactory = simIncludeGravBody.gravBodyFactory()
# setup Earth Gravity Body
earth = gravFactory.createEarth()
earth.isCentralBody = True # ensure this is the central gravitational body
mu = earth.mu
# attach gravity model to spaceCraftPlus
gravFactory.addBodiesTo(scObjectLeader)
gravFactory.addBodiesTo(scObjectFollower)
# setup MSM module
MSMmodule = msmForceTorque.MsmForceTorque()
MSMmodule.ModelTag = "msmForceTorqueTag"
scSim.AddModelToTask(dynTaskName, MSMmodule)
# define electric potentials
voltLeaderInMsgData = messaging.VoltMsgPayload()
voltLeaderInMsgData.voltage = -500 # [V] servicer potential
voltLeaderInMsg = messaging.VoltMsg().write(voltLeaderInMsgData)
voltFollowerInMsgData = messaging.VoltMsgPayload()
voltFollowerInMsgData.voltage = 500 # [V] debris potential
voltFollowerInMsg = messaging.VoltMsg().write(voltFollowerInMsgData)
# Import multi-sphere model of GOESR bus and read them into an array of strings
# For each list of 4, the first 3 values are the spacial location of an individual sphere relative to a center of
# [0,0,0] and the forth value is the radius of the sphere
path = os.path.dirname(os.path.abspath(__file__))
dataFileName = os.path.join(path, 'dataForExamples', 'GOESR_bus_80_sphs.csv')
scSphMod = open(dataFileName)
type(scSphMod)
csvreader = csv.reader(scSphMod)
rows = []
for row in csvreader:
rows.append(row)
scSphMod.close()
# Convert the strings to numbers and separate the location data from the radius data
radii = []
spherelocation = []
for row in rows:
radii.append(float(row.pop(3)))
rownum = [float(i) for i in row]
spherelocation.append(rownum)
spPosListLeader_H = spherelocation # The location of each sphere for the leader spacecraft
rListLeader = radii # radius of each sphere in the leader spacecraft
spPosListFollower_H = spherelocation # The location of each sphere for the follower spacecraft
rListFollower = radii # radius of each sphere in the follower spacecraft
# If you would like to simulate each spacecraft by a single sphere, uncomment this section (line186 - line189) of
# code and comment out the previous section lines (162-181)
# create a list of sphere body-fixed locations and associated radii using one sphere for each spacecraft
# spPosListLeader_H = [[0,0,0]] # one sphere located at origin of body frame
# rListLeader = [2] # radius of sphere is 2m
# spPosListFollower_H = [[0,0,0]] # one sphere located at origin of body frame
# rListFollower = [2] # radius of sphere is 2m
# add spacecraft to state
MSMmodule.addSpacecraftToModel(scObjectLeader.scStateOutMsg, messaging.DoubleVector(rListLeader),
unitTestSupport.npList2EigenXdVector(spPosListLeader_H))
MSMmodule.addSpacecraftToModel(scObjectFollower.scStateOutMsg, messaging.DoubleVector(rListFollower),
unitTestSupport.npList2EigenXdVector(spPosListFollower_H))
# subscribe input messages to module
MSMmodule.voltInMsgs[0].subscribeTo(voltLeaderInMsg)
MSMmodule.voltInMsgs[1].subscribeTo(voltFollowerInMsg)
# setup extForceTorque module for Leader
# the electrostatic force from the MSM module is read in through the messaging system
extFTObjectLeader = extForceTorque.ExtForceTorque()
extFTObjectLeader.ModelTag = "eForceLeader"
extFTObjectLeader.cmdForceInertialInMsg.subscribeTo(MSMmodule.eForceOutMsgs[0])
scObjectLeader.addDynamicEffector(extFTObjectLeader)
scSim.AddModelToTask(dynTaskName, extFTObjectLeader)
# setup extForceTorque module for Follower
# the electrostatic force from the MSM module is read in through the messaging system
extFTObjectFollower = extForceTorque.ExtForceTorque()
extFTObjectFollower.ModelTag = "eForceDebris"
extFTObjectFollower.cmdForceInertialInMsg.subscribeTo(MSMmodule.eForceOutMsgs[1])
scObjectFollower.addDynamicEffector(extFTObjectFollower)
scSim.AddModelToTask(dynTaskName, extFTObjectFollower)
# set initial Spacecraft States
#
# set up the Leader orbit using classical orbit elements
oe = orbitalMotion.ClassicElements()
oe.a = 42164 * 1e3 # [m] Geosynchronous Orbit
oe.e = 0.
oe.i = 0.
oe.Omega = 0.
oe.omega = 0
oe.f = 0.
r_LN, v_LN = orbitalMotion.elem2rv(mu, oe)
scObjectLeader.hub.r_CN_NInit = r_LN # m
scObjectLeader.hub.v_CN_NInit = v_LN # m/s
oe = orbitalMotion.rv2elem(mu, r_LN, v_LN)
# setup Follower object states
r_FS = np.array([0, -50.0, 0.0]) # relative position of follower, 10m behind servicer in along-track direction
r_FN = r_FS + r_LN
v_FN = v_LN
scObjectFollower.hub.r_CN_NInit = r_FN # m
scObjectFollower.hub.v_CN_NInit = v_FN # m/s
n = np.sqrt(mu / oe.a / oe.a / oe.a)
P = 2. * np.pi / n # orbit period
#
# Setup data logging before the simulation is initialized
#
numDataPoints = 1000
simulationTime = macros.sec2nano(0.1 * P)
samplingTime = simulationTime // (numDataPoints - 1)
dataRecL = scObjectLeader.scStateOutMsg.recorder()
dataRecF = scObjectFollower.scStateOutMsg.recorder()
# Add recorders to the Task
scSim.AddModelToTask(dynTaskName, dataRecL)
scSim.AddModelToTask(dynTaskName, dataRecF)
# if this scenario is to interface with the BSK Viz, uncomment the following lines
# to save the BSK data to a file, uncomment the saveFile line below
if vizSupport.vizFound:
viz = vizSupport.enableUnityVisualization(scSim, dynTaskName, [scObjectLeader, scObjectFollower]
# , saveFile=fileName,
)
#
# initialize Simulation
#
scSim.InitializeSimulation()
#
# configure a simulation stop time and execute the simulation run
#
scSim.ConfigureStopTime(simulationTime)
scSim.ExecuteSimulation()
# Retrieve the charge data of the spheres
LeaderSpCharges = unitTestSupport.columnToRowList(MSMmodule.chargeMsmOutMsgs[0].read().q)
FollowerSpCharges = unitTestSupport.columnToRowList(MSMmodule.chargeMsmOutMsgs[1].read().q)
# retrieve the logged data from the recorders
posDataL_N = dataRecL.r_BN_N
velDataL_N = dataRecL.v_BN_N
posDataF_N = dataRecF.r_BN_N
attDataL_N = dataRecL.sigma_BN
attDataF_N = dataRecF.sigma_BN
# Calculate relative position vector and magnitude in the inertial frame
relPosData_N = posDataL_N[:, 0:3] - posDataF_N[:, 0:3]
relPosMagn = np.linalg.norm(relPosData_N, axis=1)
# Project the relative position data from the inertial frame into the Hill frame of the leader spacecraft
relPosData_H = []
relXPosData_H = []
relYPosData_H = []
relZPosData_H = []
for i in range(len(relPosData_N)):
# Calculate the discrete cosine matrix for mapping from inertial frame to the Hill frame of the leader spacecraft
nrn = posDataL_N[i, :]/math.sqrt(posDataL_N[i, 0]**2 + posDataL_N[i, 1]**2 + posDataL_N[i, 2]**2)
nrh = np.cross(posDataL_N[i, 0:3], velDataL_N[i, 0:3])/np.linalg.norm(np.cross(posDataL_N[i, 0:3],
velDataL_N[i, 0:3]))
nre = np.cross(nrh, nrn)
HN = nrn, nre, nrh
# Map the relative postion data to the Hill frame of the leader spacecraft
relPosDatai_H = np.dot(HN, relPosData_N[i])
relXPosData_H.append(relPosDatai_H[0])
relYPosData_H.append(relPosDatai_H[1])
relZPosData_H.append(relPosDatai_H[2])
relPosData_H.append(relPosDatai_H)
# Collect times of each recording
timeData = dataRecL.times()
np.set_printoptions(precision=16)
figureList = plotOrbits(timeData, posDataL_N, posDataF_N, relPosMagn, attDataL_N, attDataF_N, P, spPosListLeader_H,
rListLeader, LeaderSpCharges, spPosListFollower_H, rListFollower, FollowerSpCharges,
relXPosData_H, relYPosData_H, relZPosData_H)
if show_plots:
plt.show()
# close the plots being saved off to avoid over-writing old and new figures
plt.close("all")
return figureList
def plotOrbits(timeData, posDataL_N, posDataF_N, relPosMagn, attDataL_N, attDataF_N, P, spPosListLeader_H, rListLeader,
LeaderSpCharges, spPosListFollower_H, rListFollower, FollowerSpCharges, relXPosData_H, relYPosData_H,
relZPosData_H):
#
# draw the total separation of the spacecrafts
#
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')
plt.plot(timeData * macros.NANO2SEC / P, relPosMagn[:])
plt.xlabel('Time [orbits]')
plt.ylabel('Separation [m]')
plt.title('Total separation')
figureList = {}
pltName = 'scenarioTwoChargedSC1'
figureList[pltName] = plt.figure(1)
#
# Plot relative separation in the Frame of the Leader spacecrafts
#
plt.figure(2, figsize=(5, 4))
ax = plt.axes(projection='3d')
# Set the Leader S/C as the center of the plot
r_LN_N = np.array([0., 0., 0.])
# get sphere locations
dcm_NL = RigidBodyKinematics.MRP2C(attDataF_N[0, 0:3]).transpose()
spPosL_N = np.dot(dcm_NL, np.array(spPosListLeader_H).transpose()).transpose()
radiiL = copy.deepcopy(rListLeader)
# Plot the sphere locations to model the Leader spacecraft
u = np.linspace(0, np.pi, 10)
v = np.linspace(0, 2 * np.pi, 10)
x = np.outer(np.sin(u), np.sin(v))
y = np.outer(np.sin(u), np.cos(v))
z = np.outer(np.cos(u), np.ones_like(v))
for ii in range(0, len(radiiL)):
r_SpN_N = r_LN_N + spPosL_N[ii, 0:3]
ax.plot_surface(r_SpN_N[0] + radiiL[ii] * x, r_SpN_N[1] + radiiL[ii] * y, r_SpN_N[2] + radiiL[ii] * z, color="b")
# Plot the relative position of the Follower spacecraft
ax.plot(relXPosData_H, relYPosData_H, relZPosData_H)
ax.set_xlabel('Radial(m)')
ax.set_ylabel('Along Track(m)')
ax.set_zlabel('Orbit Normal (m)')
pltName = 'scenarioTwoChargedSC2'
figureList[pltName] = plt.figure(2)
#
# Draw the sphere representation of the satellites used by the MSM in the Hill frame of the Leader spacecraft
#
SpherePlotList = SpherePlot.plotSpheres(posDataL_N, posDataF_N, attDataL_N, attDataF_N, spPosListLeader_H, rListLeader,
LeaderSpCharges, spPosListFollower_H, rListFollower, FollowerSpCharges)
figureList['scenarioTwoChargedSC3'] = SpherePlotList['Charged_Spheres']
figureList['scenarioTwoChargedSC4'] = SpherePlotList['Colorbar']
return figureList
def NormalizeData(data):
return (data - np.min(data)) / (np.max(data) - np.min(data))
#
# 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
)