Source code for test_radiation_pressure_integrated


# 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.

#
# Basilisk Integrated Test of the Solar Radiation Pressure Evaluation
#
# Purpose:  Integrated test of the spacecraft(), gravity modules and the solar
#           radiation pressure modeling.  Currently the cannonball model is only tested.
# Author:   Patrick Kenneally
# Creation Date:  June 11, 2018
#

import os

import matplotlib.pyplot as plt
import numpy as np
from Basilisk import __path__

bskPath = __path__[0]
from Basilisk.simulation import spacecraft, radiationPressure
from Basilisk.utilities import (SimulationBaseClass, macros, orbitalMotion,
                                unitTestSupport)
from Basilisk.utilities.simIncludeGravBody import gravBodyFactory


# 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() # need to update how the RW states are defined
# provide a unique test method name, starting with test_
[docs] def test_radiationPressureIntegratedTest(show_plots): """Module Unit Test""" [testResults, testMessage] = radiationPressureIntegratedTest(show_plots) assert testResults < 1, testMessage
def radiationPressureIntegratedTest(show_plots): # Create simulation variable names simTaskName = "simTask" simProcessName = "simProcess" # Create a sim module as an empty container sim = SimulationBaseClass.SimBaseClass() dynProcess = sim.CreateNewProcess(simProcessName) # create the dynamics task and specify the integration update time simulationTimeStep = macros.sec2nano(10.0) dynProcess.addTask(sim.CreateNewTask(simTaskName, simulationTimeStep)) # initialize spacecraft object and set properties scObject = spacecraft.Spacecraft() scObject.ModelTag = "spacecraftBody" sim.AddModelToTask(simTaskName, scObject) srp = radiationPressure.RadiationPressure() # default model is the SRP_CANNONBALL_MODEL srp.area = 1.0 srp.coefficientReflection = 1.3 sim.AddModelToTask(simTaskName, srp, -1) scObject.addDynamicEffector(srp) # setup Gravity Body gravFactory = gravBodyFactory() planet = gravFactory.createEarth() planet.isCentralBody = True mu = planet.mu gravFactory.createSun() spice_path = bskPath + '/supportData/EphemerisData/' gravFactory.createSpiceInterface(spice_path, '2021 MAY 04 07:47:49.965 (UTC)') gravFactory.spiceObject.zeroBase = 'Earth' sim.AddModelToTask(simTaskName, gravFactory.spiceObject, -1) srp.sunEphmInMsg.subscribeTo(gravFactory.spiceObject.planetStateOutMsgs[1]) # attach gravity model to spacecraft scObject.gravField.gravBodies = spacecraft.GravBodyVector(list(gravFactory.gravBodies.values())) # setup the orbit using classical orbit elements oe = orbitalMotion.ClassicElements() rGEO = 42000. * 1000 # meters oe.a = rGEO oe.e = 0.00001 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(mu, oe) oe = orbitalMotion.rv2elem(mu, rN, vN) # this stores consistent initial orbit elements # with circular or equatorial orbit, some angles are arbitrary print(rN) # # initialize Spacecraft States with the initialization variables # scObject.hub.r_CN_NInit = rN # m - r_BN_N scObject.hub.v_CN_NInit = vN # m/s - v_BN_N # set the simulation time n = np.sqrt(mu / oe.a / oe.a / oe.a) P = 2. * np.pi / n simulationTime = macros.sec2nano(P) # Setup data logging before the simulation is initialized numDataPoints = 100 samplingTime = simulationTime // (numDataPoints - 1) dataLog = scObject.scStateOutMsg.recorder() earthLog = gravFactory.spiceObject.planetStateOutMsgs[0].recorder() logTaskName = "logTask" dynProcess.addTask(sim.CreateNewTask(logTaskName, samplingTime)) sim.AddModelToTask(logTaskName, dataLog) sim.AddModelToTask(logTaskName, earthLog) # # initialize Simulation: This function clears the simulation log, and runs the self_init() # cross_init() and reset() routines on each module. # If the routine InitializeSimulationAndDiscover() is run instead of InitializeSimulation(), # then the all messages are auto-discovered that are shared across different BSK threads. # sim.InitializeSimulation() # # configure a simulation stop time and execute the simulation run # sim.ConfigureStopTime(simulationTime) sim.ExecuteSimulation() # unload spice kernels gravFactory.unloadSpiceKernels() # # retrieve the logged data # earthEphm = earthLog.PositionVector posData = dataLog.r_BN_N pos_rel_earth = posData[:, 0:3] - earthEphm[:, 0:3] testFailCount = 0 # zero unit test result counter testMessages = [] # create empty array to store test log messages numTruthPoints = 10 skipValue = int(len(pos_rel_earth) / (numTruthPoints - 1)) pos_rel_earth_parse = pos_rel_earth[::skipValue] # true position for un perturbed 2 body GEO orbit with cannonball SRP true_pos = np.array([[-2.18197848e+07, 3.58872415e+07, 0.00000000e+00] ,[-3.97753187e+07, 1.34888792e+07, -7.33231880e+01] ,[-3.91389859e+07, -1.52401375e+07, -3.06322198e+02] ,[-2.01838008e+07, -3.68366952e+07, -6.37764168e+02] ,[ 8.21683806e+06, -4.11950440e+07, -9.13393204e+02] ,[ 3.27532709e+07, -2.63024006e+07, -9.57828703e+02] ,[ 4.19944648e+07, 9.02522873e+05, -6.78102461e+02] ,[ 3.15828214e+07, 2.76842358e+07, -1.40473487e+02] ,[ 6.38617052e+06, 4.15047581e+07, 4.29674085e+02] ,[-2.18006914e+07, 3.58874726e+07, 7.40872311e+02]]) # compare the results to the truth values accuracy = 10.0 # meters testFailCount, testMessages = unitTestSupport.compareArray( true_pos, pos_rel_earth_parse, accuracy, "r_BN_N Vector", testFailCount, testMessages) # print out success message if no error were found if testFailCount == 0: print("PASSED ") else: print(testFailCount) print(testMessages) 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(0, 3): plt.plot(dataLog.times() * macros.NANO2SEC / P, pos_rel_earth[:, 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]') plt.title('Position Relative To Earth') if show_plots: plt.show() plt.close('all') figureList = {} fileName = os.path.basename(os.path.splitext(__file__)[0]) pltName = fileName + "srp_integrated" figureList[pltName] = plt.figure(1) return testFailCount, testMessages # # This statement below ensures that the unit test script can be run as a stand-alone python script # if __name__ == "__main__": radiationPressureIntegratedTest(True)