Source code for test_gravityGradient


# Copyright (c) 2020, 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 gravity gradient effector module.
# Author:   Hanspeter Schaub
# Creation Date:  Jan 12, 2019
#

import inspect
import os

import matplotlib.pyplot as plt
import numpy as np
import pytest
from Basilisk import __path__
from Basilisk.simulation import GravityGradientEffector
# 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 simIncludeGravBody, orbitalMotion, RigidBodyKinematics
from Basilisk.utilities import unitTestSupport

bskPath = __path__[0]

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.
[docs] @pytest.mark.parametrize("planetCase", [0, 1, 2, 3]) @pytest.mark.parametrize("cmOffset", [[[0.1], [0.15], [-0.1]], [[0.0], [0.0], [0.0]]]) # provide a unique test method name, starting with test_ def test_gravityGradientModule(show_plots, cmOffset, planetCase): r""" **Validation Test Description** This test creates a spacecraft in orbit about either Earth or Venus to check if the correct gravity gradient torque is evaluated. Multiple test scenario combinations are possible where either a single or multiple gravity bodies are included, using either zero planet ephemeris for the single planet case, or using SPICE for the multi-planet scenario. **Test Parameters** The following list discusses in detail the various test parameters used. These are test tested in all possible permutations (except show_plots of course) which is turned off for ``pytest`` usage. :param show_plots: flag to show some simulation plots :param cmOffset: center of mass offset vector in meters :param planetCase: integer flag with values (0,1,2,3). The cases consider the following simulation scenarios: - Case 0 indicates a simulation with only Earth present at (0,0,0). - Case 1 is a simulation with both Earth and Venus present using Spice, but the gravity gradient torque is only evaluated using Earth. - Case 2 is same as 1 but Venus is also included in the torque evaluation. - Case 3 is like 2 but here the spacecraft is orbiting venus. :return: None **Description of Variables Being Tested** The gravity effector torque output message is compared against a python evaluated vector. """ # each test method requires a single assert method to be called [testResults, testMessage] = run( show_plots, cmOffset, planetCase, 2.0) assert testResults < 1, testMessage
def truthGravityGradient(mu, rN, sigmaBN, hub): I = hub.IHubPntBc_B r = np.linalg.norm(rN) BN = RigidBodyKinematics.MRP2C(sigmaBN) rHatB = np.matmul(BN, rN) / r ggTorque = 3*mu/r/r/r * np.cross(rHatB, np.matmul(I, rHatB)) return ggTorque
[docs] def run(show_plots, cmOffset, planetCase, simTime): """Call this routine directly to run the unit test.""" 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)) simulationTime = macros.sec2nano(simTime) # create Earth Gravity Body gravFactory = simIncludeGravBody.gravBodyFactory() earth = gravFactory.createEarth() earth.isCentralBody = True # ensure this is the central gravitational body mu = earth.mu if planetCase: # here all planet positions are provided by Spice, and both Earth and Venus are include # for gravity acceleration calculations venus = gravFactory.createVenus() timeInitString = "2012 MAY 1 00:28:30.0" gravFactory.createSpiceInterface(bskPath + '/supportData/EphemerisData/', timeInitString, epochInMsg=True) scSim.AddModelToTask(simTaskName, gravFactory.spiceObject, -1) if planetCase == 3: # orbit should be defined relative to Venus earth.isCentralBody = False venus.isCentralBody = True mu = venus.mu gravFactory.spiceObject.zeroBase = 'venus' # spacecraft states are logged relative to Earth for plotting else: gravFactory.spiceObject.zeroBase = 'earth' # spacecraft states are logged relative to Earth for plotting # setup the orbit using classical orbit elements oe = orbitalMotion.ClassicElements() rLEO = 7000. * 1000 # meters oe.a = rLEO oe.e = 0.0001 oe.i = 33.3 * macros.D2R oe.Omega = 48.2 * macros.D2R oe.omega = 347.8 * macros.D2R oe.f = 85.3 * 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 # setup basic spacecraft module scObject = spacecraft.Spacecraft() scObject.ModelTag = "bskTestSat" IIC = [[500., 0., 0.] , [0., 800., 0.] , [0., 0., 350.]] scObject.hub.r_BcB_B = cmOffset scObject.hub.mHub = 100.0 # kg - spacecraft mass scObject.hub.IHubPntBc_B = IIC scObject.hub.r_CN_NInit = rN # m - r_BN_N scObject.hub.v_CN_NInit = vN # m/s - v_BN_N scObject.hub.sigma_BNInit = [[0.1], [0.2], [-0.3]] # sigma_BN_B scObject.hub.omega_BN_BInit = [[0.0], [0.0], [0.0]] # rad/s - omega_BN_B scSim.AddModelToTask(simTaskName, scObject) scObject.gravField.gravBodies = spacecraft.GravBodyVector(list(gravFactory.gravBodies.values())) # add gravity gradient effector ggEff = GravityGradientEffector.GravityGradientEffector() ggEff.ModelTag = scObject.ModelTag ggEff.addPlanetName(earth.planetName) if planetCase >= 2: ggEff.addPlanetName(venus.planetName) scObject.addDynamicEffector(ggEff) scSim.AddModelToTask(simTaskName, ggEff) # # Setup data logging before the simulation is initialized # numDataPoints = 50 samplingTime = unitTestSupport.samplingTime(simulationTime, simulationTimeStep, numDataPoints) dataLog = scObject.scStateOutMsg.recorder(samplingTime) dataLogGG = ggEff.gravityGradientOutMsg.recorder(samplingTime) scSim.AddModelToTask(simTaskName, dataLog) scSim.AddModelToTask(simTaskName, dataLogGG) # # 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 attData = dataLog.sigma_BN ggData = dataLogGG.gravityGradientTorque_B np.set_printoptions(precision=16) # # plot the results # if show_plots: plt.close("all") # clears out plots from earlier test runs # draw the inertial position vector components plt.close("all") # clears out plots from earlier test runs plt.figure(1) for idx in range(0, 3): plt.plot(dataLog.times() * macros.NANO2MIN, attData[:, idx], color=unitTestSupport.getLineColor(idx, 3), label=r'$\sigma_' + str(idx) + '$') plt.legend(loc='lower right') plt.xlabel('Time [min]') plt.ylabel(r'MRP Attitude $\sigma_{B/N}$') plt.figure(2) for idx in range(0, 3): plt.plot(dataLog.times() * macros.NANO2MIN, posData[:, idx]/1000, color=unitTestSupport.getLineColor(idx, 3), label=r'$r_' + str(idx) + '$') plt.legend(loc='lower right') plt.xlabel('Time [min]') plt.ylabel(r'Inertial Position coordinates [km]') plt.figure(3) for idx in range(0, 3): plt.plot(dataLogGG.times() * macros.NANO2MIN, ggData[:, idx] , color=unitTestSupport.getLineColor(idx, 3), label=r'$r_' + str(idx) + '$') plt.legend(loc='lower right') plt.xlabel('Time [min]') plt.ylabel(r'GG Torque [Nm]') plt.show() plt.close("all") # compare gravity gradient torque vector to the truth accuracy = 1e-10 for rV, sV, ggV in zip(posData, attData, ggData): ggTruth = truthGravityGradient(mu, rV[0:3], sV[0:3], scObject.hub) testFailCount, testMessages = unitTestSupport.compareVector(ggV[0:3], ggTruth, accuracy, "gravityGradientTorque_B", testFailCount, testMessages) print("Accuracy used: " + str(accuracy)) if testFailCount == 0: print("PASSED: Gravity Effector") else: print("Failed: Gravity Effector") return testFailCount, testMessages
# close the plots being saved off to avoid over-writing old and new figures if __name__ == '__main__': run(True, # show_plots [[0.0], [0.0], [0.0]], # cmOffset 3, # planetCase 3600) # simTime (seconds)