Source code for scenarioSatelliteConstellation

#
#  ISC License
#
#  Copyright (c) 2023, 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 scenario demonstrates how to set up a Walker-Delta constellation of satellites. Note that this scenario 
uses the stand-alone Basilisk architecture rather than using the ''examples/FormationBskSim`` or ``examples/MultiSatBskSim``
architectures for simultaneously simulating multiple spacecraft.

The script is found in the folder ``basilisk/examples`` and executed by using::

    python3 scenarioSatelliteConstellation.py

When designing a satellite constellation, symmetry provides repeatable performance and makes design and operations simpler. One 
such symmetric constellation design methodology is the Walker constellation which consists of circular orbits in evenly spaced 
planes at a chosen inclination. Walker constellation layouts are fully specified by four parameters written as "i:T/P/F" where:

    * i = inclination angle
    * T = total number of satellites in the constellation
    * P = numer of orbit planes evenly spaced in node angle
    * F = relative spacing between satellites in adjacent orbit planes

In order to simulate the constellation a semi-major axis 'a' is also specified as an input. The total number of satellites T must 
be specified as an integer multiple of the number of planes P. The relative spacing F must be an integer between 0 and (P-1).

Illustration of Simulation Results
----------------------------------

::

    show_plots = True, a = 29994000, i = 56, T = 24, P = 3, F = 1

The default scenario is modelled after the Galileo constellation with Walker numbers 56:24/3/1 and a semi-major axis of 29,994 km.

.. image:: /_images/Scenarios/scenarioSatelliteConstellation.svg
   :align: center

Each line color corresponds to a different orbital plane, in this case 3 planes. Green circles mark the start of each of the 24 
individual satellite's ground tracks and a red circle marks their end.

Simulation Visualization In Vizard
----------------------------------

To add a visualization of each spacecraft's Earth nadir pointing, a cone aligned with the pointing direction
is created through the ``vizInterface``::

    vizSupport.createConeInOut(viz, fromBodyName=scList[i].ModelTag, toBodyName='earth', coneColor='teal',
                               normalVector_B=[1, 0, 0], incidenceAngle=30/macros.D2R, isKeepIn=False,
                               coneHeight=10.0, coneName='pointingCone')

.. image:: /_images/static/scenarioSatelliteConstellationVizard.jpg
   :scale: 50%
   :align: center

"""

#
# Basilisk Scenario Script and Integrated Test
#
# Purpose:  Demonstrates how to use the stand alone Basilisk architecture to automate 
#           the setup of a Walker constellation.
# Author:   Andrew Morell
# Creation Date:  Dec. 13, 2023
#

import os

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

bskPath = __path__[0]
fileName = os.path.basename(os.path.splitext(__file__)[0])

from Basilisk.simulation import (spacecraft, ephemerisConverter)
from Basilisk.utilities import (SimulationBaseClass, macros, orbitalMotion,
                                simIncludeGravBody, unitTestSupport, vizSupport)
from Basilisk.fswAlgorithms import locationPointing
from Basilisk.simulation import (simpleNav, planetEphemeris)

# always import the Basilisk messaging support

[docs] def run(show_plots, a, i, T, P, F): """ The scenario can be run with the followings setups parameters: Args: show_plots (bool): Determines if the script should display plots a (float): semi-major axis of all satellites [m] i (float): orbital inclination [rad] T (int): total number of satellites P (int): number of equally spaced orbital planes F (int): phasing between satellites in adjacent places """ # Create simulation variable names simTaskName = "simTask" simProcessName = "simProcess" # Create a sim module as an empty container scSim = SimulationBaseClass.SimBaseClass() # (Optional) If you want to see a simulation progress bar in the terminal window, the # use the following SetProgressBar(True) statement scSim.SetProgressBar(True) # create the simulation process dynProcess = scSim.CreateNewProcess(simProcessName) # create the dynamics task and specify the integration update time simulationTimeStep = macros.sec2nano(30.) dynProcess.addTask(scSim.CreateNewTask(simTaskName, simulationTimeStep)) # 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 # Create the ephemeris converter module spiceObject = gravFactory.createSpiceInterface(time="2029 June 12 5:30:30.0", epochInMsg=True) spiceObject.zeroBase = 'Earth' scSim.AddModelToTask(simTaskName, spiceObject) ephemObject = ephemerisConverter.EphemerisConverter() ephemObject.ModelTag = "ephemData" ephemObject.addSpiceInputMsg(spiceObject.planetStateOutMsgs[0]) # Earth scSim.AddModelToTask(simTaskName, ephemObject) # spacecraft inertia I = [900., 0., 0., 0., 800., 0., 0., 0., 600.] # set the simulation time to be two orbit periods n = np.sqrt(mu / a / a / a) Pr = 2. * np.pi / n simulationTime = macros.sec2nano(1.5 * Pr) # create a logging task object of the spacecraft output message at the desired down sampling ratio numDataPoints = 500 samplingTime = unitTestSupport.samplingTime(simulationTime, simulationTimeStep, numDataPoints) if (np.mod(T,1) !=0 or np.mod(P,1) !=0 or np.mod(F,1) != 0 or T < 1 or P < 1 or F < 1): raise Exception('number of satellites T, number of planes P, and relative spaceing F must be positive integer numbers') if (np.mod(T,P) !=0): raise Exception('number of satellites T must be an integer multiple of P') # pre-compute Walker constellation parameters PU = 360/T # patter unit in degrees S = T/P # number of satellites in each plane delta_RAAN = S*PU # delta RAAN in degrees delta_nu = P*PU # in plane spacing of satellites oe = orbitalMotion.ClassicElements() oe.a = a oe.e = 0.01 oe.i = i # initialize T spacecraft objects and set properties scList = [] navList = [] guideList = [] dataRec = [] for i in range(T): # initialize spacecraft object scList.append(spacecraft.Spacecraft()) scList[i].ModelTag = "bsk-Sat-"+str(i) scSim.AddModelToTask(simTaskName, scList[i]) # attach gravity model to spacecraft scList[i].gravField.gravBodies = spacecraft.GravBodyVector(list(gravFactory.gravBodies.values())) # setup the orbit using classical orbit elements RAAN = delta_RAAN * np.floor(i/S) oe.Omega = RAAN * macros.D2R oe.omega = F * np.floor(i/S) * macros.D2R oe.f = delta_nu * np.mod(i,S) * macros.D2R rN, vN = orbitalMotion.elem2rv(mu, oe) scList[i].hub.r_CN_NInit = rN # [m] scList[i].hub.v_CN_NInit = vN # [m/s] # set the spacecraft mass and inertia scList[i].hub.mHub = 750.0 # kg - spacecraft mass scList[i].hub.IHubPntBc_B = unitTestSupport.np2EigenMatrix3d(I) # set up prescribed attitude guidance navList.append(simpleNav.SimpleNav()) navList[i].ModelTag = "SimpleNav"+str(i) scSim.AddModelToTask(simTaskName, navList[i]) navList[i].scStateInMsg.subscribeTo(scList[i].scStateOutMsg) # setup nadir pointing guidance module guideList.append(locationPointing.locationPointing()) guideList[i].ModelTag = "nadirPoint"+str(i) guideList[i].pHat_B = [1, 0, 0] guideList[i].scTransInMsg.subscribeTo(navList[i].transOutMsg) guideList[i].scAttInMsg.subscribeTo(navList[i].attOutMsg) guideList[i].celBodyInMsg.subscribeTo(ephemObject.ephemOutMsgs[0]) scSim.AddModelToTask(simTaskName, guideList[i]) scList[i].attRefInMsg.subscribeTo(guideList[i].attRefOutMsg) # record spacecraft states dataRec.append(scList[i].scStateOutMsg.recorder(samplingTime)) scSim.AddModelToTask(simTaskName, dataRec[i]) # If you wish to transmit the simulation data to the United based Vizard Visualization application, # then uncomment the following saveFile line if vizSupport.vizFound: viz = vizSupport.enableUnityVisualization(scSim, simTaskName, scList, # saveFile=__file__ ) for i in range(T): vizSupport.createConeInOut(viz, fromBodyName=scList[i].ModelTag, toBodyName='earth', coneColor='teal', normalVector_B=[1, 0, 0], incidenceAngle=30/macros.D2R, isKeepIn=False, coneHeight=10.0, coneName='pointingCone') # initialize simulation scSim.InitializeSimulation() # configure a simulation stop time and execute the simulation run scSim.ConfigureStopTime(simulationTime) scSim.ExecuteSimulation() # plot the ground track plt.close("all") # clears out plots from earlier test runs figureList = {} plt.figure(1) for n in range(T): # calculate latitude and longitude [lon, lat] = rv2latlon(dataRec[n]) color = unitTestSupport.getLineColor(int(np.floor(n/S)), 3) plt.scatter(lat, lon, s=5, c=[color]) plt.scatter(lat[0], lon[0], s=20, c='#00A300') plt.scatter(lat[-1], lon[-1], s=20, c='#A30000') plt.xlim([-180,180]) plt.ylim([-90,90]) plt.xticks(range(-180,181,30)) plt.yticks(range(-90,91,15)) plt.xlabel('longitude (degrees)') plt.ylabel('lattitude (degrees)') pltName = fileName figureList[pltName] = plt.figure(1) 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 rv2latlon(dataRec): times = dataRec.times() * macros.NANO2SEC rECI = dataRec.r_BN_N lat = np.zeros(times.shape) lon = np.zeros(times.shape) for i in range(len(times)): theta = times[i]*planetEphemeris.OMEGA_EARTH dcm_ECI2ECEF = [[np.cos(theta), np.sin(theta), 0], [-np.sin(theta), np.cos(theta), 0], [0, 0, 1]] rECEF = dcm_ECI2ECEF@rECI[i,:] lat[i] = np.arcsin(rECEF[2]/np.linalg.norm(rECEF)) * macros.R2D lon[i] = np.arctan2(rECEF[1],rECEF[0]) * macros.R2D return lat, lon if __name__ == "__main__": run( True, # show_plots 29994000, # semi-major axis [m] 56, # orbit inclination [deg] 24, # total number of satellites (int) 3, # number of orbit planes (int) 1 # phasing (int) )