High Fidelity Orbital Dynamics

[AmolikaSoni]

Hi @schaubh ,
I have been using the SC Object Position & Velocity to mimic GPS data and subsequently use it to test my orbit determination algorithm for a SSO orbit small sat. I tried to compare the dynamics with the library offered by https://www.orekit.org/

I observed that the orbit data generated by https://www.orekit.org/ and BSK diverge and reach up to an error of ~2kms over 24hr of orbit.
I have plots to demonstrate the same.

I was wondering if the orbital dynamics in Basilisk can be updated to reach this level of fidelity. I see that the gravity data being used it GGM03S https://www2.csr.utexas.edu/grace/gravity/ggm03/. Upon a bit of Google search, I found that GGM05 is also available now https://www2.csr.utexas.edu/grace/gravity/

Are there any plans to update the gravity data?
Is there any provision to customize and use higher-fidelity data?

[AmolikaSoni]

[bbercoviciUspace]

Hello @AmolikaSoni ,

There are a number of things that may have come into play to explain the divergence.

Have you taken steps to ensure that the BSK scenario you ran effectively makes use of the full spherical harmonics model (and doesn't default to point-mass gravity) ? Something else to check would be the proper modeling of the Earth's rotation : if I am not mistaken, unless the Earth gravity model is explicitly connected to Spice, the Earth will basically remain inertially fixed.

I would also try to display the osculating orbital elements instead of the cartesian components of position/velocity to see if something obvious stands out.

[schaubh]

Howdy @AmolikaSoni , I agree with @bbercoviciUspace comments. Basilisk dynamics includes the ability to include N-th order spherical harmonics, polyhedral gravity models, lumped mass models. By default, if you just use createEarth you get a point mass model with a non-rotating object. To give the planet a specific orientation that varies with time, you can connect to the Spice interface as discussed in the examples.

You ask if we are going to include higher fidelity gravity modeling, we already do this. We just don't provide the data files for higher order spherical harmonics. You can readily create this on your own from the sources you site. Take a look at

https://avslab.github.io/basilisk/_modules/scenarioOrbitMultiBody.html

to see how a spherical harmonics file is loaded.

[AmolikaSoni]

Hi @schaubh & @bbercoviciUspace ,
Thanks for your response. I believe I am doing all the things that you have mentioned.

• Creating sim base class
• Create SC object and assign highest priority
• Create GravBodyFactory
• Create exponential atmosphere, drag effector, eclipse, and solar flux objects (just to ensure all the forces from the environment are accounted for)

Here is a snippet of my python code. Can you please take a look and let me know if I have missed anything?

    # 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(1)
    dynProcess.addTask(scSim.CreateNewTask(simTaskName, simulationTimeStep))

    # setup the simulation tasks/objects
    # initialize spacecraft object and set properties
    # The dynamics simulation is setup using a Spacecraft() module.
    scObject = spacecraft.Spacecraft()
    scObject.ModelTag = "bsk-Sat"

    

    # setup Gravity Body
    gravFactory = simIncludeGravBody.gravBodyFactory()
    gravBodies = gravFactory.createBodies('earth', 'mars barycenter', 'sun', 'moon', "jupiter barycenter")
    gravBodies['earth'].isCentralBody = True

    gravBodies['earth'].useSphericalHarmonicsGravityModel(bskPath + '/supportData/LocalGravData/GGM03S.txt', 100)

    mu = gravBodies['earth'].mu
    r_eq = gravBodies['earth'].radEquator
    sun = 0
    earth = 1

    gravFactory.addBodiesTo(scObject)
    scObject.gravField.gravBodies = spacecraft.GravBodyVector(list(gravFactory.gravBodies.values()))
    timeInitString = "2025 JAN 15 00:00:00.000"
    spiceObject = gravFactory.createSpiceInterface(time=timeInitString, epochInMsg=True)
    spiceObject.zeroBase = 'Earth'



    gravFactory.spiceObject.zeroBase = 'Earth'
    EarthEphemObject = ephemerisConverter.EphemerisConverter()
    EarthEphemObject.addSpiceInputMsg(gravFactory.spiceObject.planetStateOutMsgs[sun])
    EarthEphemObject.addSpiceInputMsg(gravFactory.spiceObject.planetStateOutMsgs[earth])

    pyswice.furnsh_c(gravFactory.spiceObject.SPICEDataPath + 'de430.bsp')  # solar system bodies
    pyswice.furnsh_c(gravFactory.spiceObject.SPICEDataPath + 'naif0012.tls')  # leap second file
    pyswice.furnsh_c(gravFactory.spiceObject.SPICEDataPath + 'de-403-masses.tpc')  # solar system masses
    pyswice.furnsh_c(gravFactory.spiceObject.SPICEDataPath + 'pck00010.tpc')  # generic Planetary Constants
    
    #Atmosphere
    expAtmo = exponentialAtmosphere.ExponentialAtmosphere()
    expAtmo.ModelTag = "ExpAtmo"
    expAtmo.addSpacecraftToModel(scObject.scStateOutMsg)
    expAtmo.planetRadius = r_eq
    # Drag
    dragEff  = dragDynamicEffector.DragDynamicEffector()
    dragEff.ModelTag = "DragEff"
    dragEff.coreParams.projectedArea = .5
    dragEff.coreParams.dragCoeff = 2.2
    dragEff.coreParams.comOffset = [0., 0., 0.]
    scObject.addDynamicEffector(dragEff)
    dragEff.atmoDensInMsg.subscribeTo(expAtmo.envOutMsgs[0])
    #Eclipse
    eclipseObject  = eclipse.Eclipse()
    """
    Specify what celestial object is causing an eclipse message.
    """
    eclipseObject.ModelTag = "eclipseObject"
    eclipseObject.sunInMsg.subscribeTo(gravFactory.spiceObject.planetStateOutMsgs[sun])
    # add all celestial objects in spiceObjects except for the sun (0th object)
    for c in range(1, len(gravFactory.spiceObject.planetStateOutMsgs)):
        eclipseObject.addPlanetToModel(gravFactory.spiceObject.planetStateOutMsgs[c])
    eclipseObject.addSpacecraftToModel(scObject.scStateOutMsg)
    
    #Solar Flux
    sf = solarFlux.SolarFlux()
    sf.sunPositionInMsg.subscribeTo(gravFactory.spiceObject.planetStateOutMsgs[sun])
    sf.spacecraftStateInMsg.subscribeTo(scObject.scStateOutMsg)
    sf.eclipseInMsg.subscribeTo(eclipseObject.eclipseOutMsgs[0])
    sf.ModelTag = "SolarFlux"
    
    # add spacecraft object to the simulation process
    scSim.AddModelToTask(simTaskName, scObject,1000)
    scSim.AddModelToTask(simTaskName, spiceObject, 999)
    scSim.AddModelToTask(simTaskName, gravFactory.spiceObject, None, 998)
    scSim.AddModelToTask(simTaskName, expAtmo, None, 996)
    scSim.AddModelToTask(simTaskName, dragEff, None, 995)
    scSim.AddModelToTask(simTaskName, eclipseObject, None, 994)
    scSim.AddModelToTask(simTaskName, sf, None, 992)
    scSim.AddModelToTask(simTaskName, EarthEphemObject, 989)

    #
    #   setup orbit 
    #

    # To set the spacecraft initial conditions, the following initial position and velocity variables are set:
    rN = np.array([-267322.99915222195,6962999.124678645,8.149072527885437e-10])
    vN = np.array([1019.2392593879555,31.101414308892743,7494.238898865543])
    scObject.hub.r_CN_NInit = rN  # m   - r_BN_N
    scObject.hub.v_CN_NInit = vN  # m/s - v_BN_N
    oe = orbitalMotion.rv2elem(mu, rN, vN)      # this stores consistent initial orbit elements

[bbercoviciUspace]

I think you need to add the following

gravBodies['earth'].useSphericalHarmonicsGravityModel = True

My guess is that that above scenario only accounts for a point-mass earth due to the missing line. My recommendation to plot the orbital elements still holds : I think you will find that the RAAN is constant, whereas the high-fidelity data you are comparing the BSK orbit with must exhibit a secular drift caused by J2.

[AmolikaSoni]

@bbercoviciUspace I added the line as you have mentioned. (Although, I did not find any such line in the example codes and hence I did not add it to my original code).
Nevertheless, here is a graph of all the orbit parameters.

[bbercoviciUspace]

How do these orbital elements compare to the reference solution ? Do they look the same when useSphericalHarmonicsGravityModel is left unset ?

[AmolikaSoni]

Hi @schaubh, I ran the scenarioOrbitMultiBody.py example as well for a longer duration. The errors between spice and basilisk position go up to ~15kms. Is that expected?

[AmolikaSoni]

@bbercoviciUspace I did the comparison for orbital elements as you suggested. There is absolutely no change.

[schaubh]

Howdy @AmolikaSoni , you create an exponential atmosphere module, but you never specify the atmosphere parameters? You only set req. Given this is a LEO satellite, this could explain the difference. The other values are defaulted to zero.

There is basilisk/utilities/simSetPlanetEnvironment.py helper functions that will set default exponential atmosphere values for earth. The module PDF documentation talks about this. see https://avslab.github.io/basilisk/Documentation/utilities/simSetPlanetEnvironment.html.

Right now we just have sample values for Earth. If you are comparing BSK to another simulation, I'd ensure you are using the same atmosphere model and values.

In your sim you are making an eclipse module, but it is not used anywhere. The solar flux module just provides the local solar value in the solar system. Given you are doing a test at LEO, I think the SRP force would be comparatively small. Right now you are simply not setting up your atmosphere model.