GMAT Code Example

danger

This document contains confidential data and is not intended to be shared outside of SAGE!

Authors	Julien Purificato
Responsible	Miles Timpe
Last Updated	10/12/2025, 12:48:54 PM
Last Author	Kai Berszin

What's the purpose of this code

This code was developed to provide an example of how GMAT can be used to simulate both an orbit (the MOPS team's main task) and attitude (the ADCS team's main task). The task is to simulate a satellite's orbit and orient it toward the ground station when there is access, and toward the Sun otherwise. The code can also be used to compare orbit propagation results as predicted by GMAT with the predictions from other software such as OSTK. It can be found in sage-flight-dynamics\flight_dynamics\orbit_propagation\gmat under the name DetectAndPointGMAT.script or in the github: https://github.com/aris-space/sage-flight-dynamics

The code primarily takes three types of input:

• Orbit characteristics (altitude, local time at descending node, etc.) that you can change within the code after the creation of the Spacecraft GEOSat.
• Simulation times (initial date-time, timestep, end instant) that you can change within the code
• Information on the ground station (name, latitude, longitude, altitude, etc.) that can be changed after the creation of the GroundStation myStation.

In return, it creates one CSV file containing all the valuable datas OSTK was already giving:

• Timestamps (UTC). ✔️ The timestamps are given in TAI which is UTC + 37s, thus the initial timestep used to compare the initial datas are displayed as 12:00:00 in the OSTK CSV and 12:00:37 in the GMAT CSV
• Julian Date (UTC). The UTCModJulian of GMAT (e.g 30683.0947573436) corresponds to the MJD of OSTK (e.g 60682.5947569445) with a constant offset of -29999.5.✔️
• 3D Positions (GCRF = Earth-centered, non-rotating inertial frame) ✔️
• 3D Velocities (GCRF) ✔️
• Latitude, Longitude, and Altitude ✔️
• During accesses: Azimuth, Elevation, and Range ✔️
• Yaw, Pitch, and Roll : they are rotations around the z-, y-, and x-axes respectively, in the VVLH frame (x is orbital velocity, z is nadir) ? The nadir pointing seems to correctly orient the z-axis of the satellite, but not the x or y, which makes it impossible to check if this step is done correctly (dynamically point to Earth and observe zero pitch and roll in a VVLH frame). I kept in the script the definition of the VVLH frame, an example of how to do Sun-pointing, and the reporting of the three Euler angles and the quaternions.
• Indication on whether the satellite is in eclipse ✔️ and/or in access ✔️. While OSTK provides a boolean, GMAT creates report files ("EclipseLocator1.txt" and "ContactLocator1.txt" respectively) which indicate timestamps of eclipse and access.
• During accesses: Acquisition of Signal (AOS) ✔️, Time of Closest Approach (TCA) ❌ , Loss of Signal (LOS), and duration of the access ✔️

Discussion of results

Comparison of the evolution of orbit between GMAT & OSTK is given below: /TBD/

The inital orbital elements created in OSTK can be exactly reproduced in GMAT using the following code:

import os
import math
import numpy as np
import pandas as pd

import plotly.express as px
import plotly.express as px
import pandas as pd
import numpy as np

from ostk.core.filesystem import Directory

from ostk.mathematics.geometry.d3.object import Cuboid
from ostk.mathematics.geometry.d3.object import Composite
from ostk.mathematics.geometry.d3.object import Point

from ostk.physics import Environment
from ostk.physics.coordinate import Frame
from ostk.physics.environment.atmospheric import Earth as EarthAtmosphericModel
from ostk.physics.environment.gravitational import Earth as EarthGravitationalModel
from ostk.physics.environment.magnetic import Earth as EarthMagneticModel
from ostk.physics.environment.object.celestial import Earth
from ostk.physics.environment.object.celestial import Moon
from ostk.physics.environment.object.celestial import Sun
from ostk.physics.time import DateTime
from ostk.physics.time import Duration
from ostk.physics.time import Instant
from ostk.physics.time import Scale
from ostk.physics.time import Time
from ostk.physics.unit import Length
from ostk.physics.unit import Mass

from ostk.astrodynamics import Dynamics
from ostk.astrodynamics.dynamics import Thruster
from ostk.astrodynamics.event_condition import RealCondition
from ostk.astrodynamics.event_condition import LogicalCondition
from ostk.astrodynamics.event_condition import BooleanCondition
from ostk.astrodynamics.trajectory.orbit.model.brouwerLyddaneMean import (
    BrouwerLyddaneMeanShort,
)
from ostk.astrodynamics.trajectory import Orbit
from ostk.astrodynamics.trajectory import Sequence
from ostk.astrodynamics.trajectory import StateBuilder
from ostk.astrodynamics.trajectory.state import NumericalSolver
from ostk.astrodynamics.trajectory.state import CoordinateBroker, CoordinateSubset
from ostk.astrodynamics.trajectory.state.coordinate_subset import CartesianPosition
from ostk.astrodynamics.trajectory.state.coordinate_subset import CartesianVelocity
from ostk.astrodynamics.trajectory.orbit.model.kepler import COE
from ostk.astrodynamics.flight.system import SatelliteSystem
from ostk.astrodynamics.flight.system import PropulsionSystem
from ostk.astrodynamics.guidance_law import QLaw
from ostk.mathematics.object import RealInterval
from ostk.mathematics.geometry import Angle as MathAngle
from ostk.mathematics.geometry.d3.transformation.rotation import Quaternion
from ostk.mathematics.geometry.d3.transformation.rotation import EulerAngle

from ostk.physics import Environment
from ostk.physics.unit import Length
from ostk.physics.unit import Angle
from ostk.physics.time import Scale
from ostk.physics.time import Instant
from ostk.physics.time import Duration
from ostk.physics.time import Interval
from ostk.physics.time import DateTime
from ostk.physics.time import Time
from ostk.physics.coordinate import Frame
from ostk.physics.coordinate import Transform
from ostk.physics.coordinate.spherical import LLA
from ostk.physics.coordinate.frame.provider import Dynamic as DynamicFrameProvider

from ostk.astrodynamics import Trajectory
from ostk.astrodynamics.trajectory import Orbit
from ostk.astrodynamics.flight import Profile
from ostk.astrodynamics.flight.profile.model import Transform as FlightProfileTransform
from ostk.astrodynamics.access import Generator as AccessGenerator
from ostk.astrodynamics.access import AccessTarget
from ostk.astrodynamics.access import VisibilityCriterion
from ostk.astrodynamics.viewer import Viewer
from ostk.astrodynamics.trajectory.orbit.model import SGP4
from ostk.astrodynamics.trajectory.orbit.model.sgp4 import TLE
from ostk.astrodynamics.utilities import convert_state
from ostk.astrodynamics.utilities import compute_trajectory_geometry
from ostk.astrodynamics.utilities import compute_time_lla_aer_coordinates

#1) create the orbit path and find the accesses in a time range (some modification from orbit_and_accesses_example.py)

# User-defined functions
def compute_access_geometry(access): #compute AER (Azimuth-Elevation-Range) and LLA (Latitude-Longitude-Altitude) for an entire access
    return [
        compute_time_lla_aer_coordinates(state, access_target.get_position(), environment)
        for state in satellite_orbit.get_states_at(
            access.get_interval().generate_grid(Duration.seconds(1.0))
        )
    ]

def fix_timestamp(s):
    """Fix timestamp formatting bug currently present in compute_accesses() function."""    
    s = str(s)
    return s.split('.')[0]+'.'+''.join(s.split('.')[1:])


def in_eclipse(state): #Checks whether the satellite is in Eclipse at a given states (position,velocity,timestamps)
    # Global environment
    environment.set_instant(state.get_instant())
    return environment.is_position_in_eclipse(state.get_position())


if __name__=="__main__":

    ### Define simulation environment

    # Set simulation environment to include Earth and the Sun
    environment = Environment.default()
    earth = environment.access_celestial_object_with_name("Earth")
    sun = environment.access_celestial_object_with_name("Sun")


    # Set initial epoch
    epoch = Instant.date_time(DateTime(2025, 1, 7, 12, 0, 0), Scale.UTC)

    ### Satellite orbit

    # Satellite IDs can be used to generalize this script to multiple satellites
    sat_id = "SAGE CubeSat"

    # Define a simple sun-synchronous orbit (SSO) for SAGE
    satellite_orbit = Orbit.sun_synchronous(
        epoch=epoch,
        altitude=Length.kilometers(510.0),
        local_time_at_descending_node=Time(18, 0, 0),
        celestial_object=earth,
        argument_of_latitude=Angle.degrees(90.0)
    )
kepler_model = satellite_orbit.access_kepler_model()
keplerian_elements = kepler_model.get_classical_orbital_elements()
print(keplerian_elements)

#If the print doesn't work
import sys
sys.stdout.flush()

print(keplerian_elements)

and gives the following elements:

With this orbit we observe that 3D positions, 3D Velocities and LLA gives the exact same values (less than 1e-3 % of relative difference) at initial time. The Azimuth, Elevation and range are calculated withing the code as these variables are not implemented in GMAT. The computation method was verified using matlab and the differences in the values between OSTK and GMAT arises from slightly (/Quantify/ /TBD/) different values in the position throughout time simulation.

The pointing part need running the code 2 times : first to find the access epoch and second request change of pointing during the obtained epochs.

What should be improved and what still need to be done

• Further comparison of values between the two software
• Check correct implementation of YPR data (I have also tried forcing the quaternion but that doesn't work either)
• Find a way to only run the code once
• Make Station Pointing Changelog

Date	Revision	Change
2025/07/15	01	Initial commit
2025/08/01	02	Final commit