Add Walker Delta constellation setup scenario

docs/source/Support/bskReleaseNotes.rst

@@ -32,6 +32,7 @@ Basilisk Release Notes
 
 Version |release|
 -----------------
+- Created a new example scenario :ref:`scenarioSatelliteConstellation` demonstrating setup of a Walker-Delta constellation
 - Created a new :ref:`pinholeCamera` module to support generation of landmarks-based measurements around a
   small body.
 - Corrected a memory leak in the ``swig`` access to standard vectors inside messages.

examples/_default.rst

@@ -294,6 +294,7 @@ Stand Alone Architecture
    Electrostatic Tractor Debris Reorbiting <scenarioDebrisReorbitET>
    Attitude-Driven differential drag control <scenarioDragRendezvous>
    Servicer approaching a debris object with 3 flight modes <scenarioRendezVous>
+   Walker-Delta Satellite Constellation <scenarioSatelliteConstellation>
 
 
 ``FormationBskSim`` Architecture

examples/scenarioSatelliteConstellation.py

@@ -0,0 +1,292 @@
+#
+#  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 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
+
+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 an altitude of 29,994 km, inclination of 
+56 degrees, and 24 satellites across 3 planes with a phasing offset of 1.
+
+.. image:: /_images/Scenarios/scenarioSatelliteConstellationGroundTrack.svg
+   :align: center
+
+Simulation Visualization In Vizard
+----------------------------------
+
+To add a visualization of each spacecraft's Earth nadir pointing, a cone algined 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/Scenarios/scenarioSatelliteConstellationVizard.svg
+   :align: center
+
+"""
+
+#
+# Basilisk Scenario Script and Integrated Test
+#
+# Purpose:  Integrated test of the spacecraft() and gravity modules.  Illustrates
+#           a 3-DOV spacecraft on a range of orbit types.
+# Author:   Andrew Morell
+# Creation Date:  Dec. 13, 2023
+#
+
+import os
+from copy import copy
+
+import matplotlib.pyplot as plt
+import numpy as np
+# To play with any scenario scripts as tutorials, you should make a copy of them into a custom folder
+# outside of the Basilisk directory.
+#
+# To copy them, first find the location of the Basilisk installation.
+# After installing, you can find the installed location of Basilisk by opening a python interpreter and
+# running the commands:
+from Basilisk import __path__
+
+bskPath = __path__[0]
+fileName = os.path.basename(os.path.splitext(__file__)[0])
+
+# Copy the folder `{basiliskPath}/examples` into a new folder in a different directory.
+# Now, when you want to use a tutorial, navigate inside that folder, and edit and execute the *copied* integrated tests.
+
+
+# import simulation related support
+from Basilisk.simulation import (spacecraft, ephemerisConverter)
+# general support file with common unit test functions
+# import general simulation support files
+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
+
+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
+        i (float): orbital inclination
+        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(10.)
+    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 line
+    viz = vizSupport.enableUnityVisualization(scSim, simTaskName, scList,
+                                              saveFile=__file__
+                                              # liveStream=True
+                                              )
+    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
+    for n in range(T):
+        # calculate lattitude 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='#750000')
+        plt.scatter(lat[-1], lon[-1], s=20, c='#00A300')
+    plt.xlim([-180,180])
+    plt.ylim([-90,90])
+    plt.xticks(range(-180,200,20))
+    plt.yticks(range(-90,105,15))
+    plt.xlabel('longitude (degrees)')
+    plt.ylabel('latitude (degrees)')
+    img = plt.imread("docs/source/_images/Scenarios/scenarioSatConstellationMercator.png")
+    plt.imshow(img, extent=[-180,180,-90,90])
+
+    if show_plots:
+        plt.show()
+
+    # close the plots being saved off to avoid over-writing old and new figures
+    plt.close("all")
+
+    return
+
+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)
+    )