1. Unperturbed Earth-orbiting Satellite

The example described on this page is that of Asterix, a single satellite that is orbiting the Earth. The code for this tutorial is given here on Github, and is also located in your tudat bundle at:

tudatBundle/tudatExampleApplications/satellitePropagatorExamples/SatellitePropagatorExamples/singleSatellitePropagator.cpp

For this example, we have the following problem statement:

Given the position and velocity of the Asterix satellite at a certain point in time with respect to the Earth, what will its position and velocity be after a Julian day has passed?

There are a number of assumptions we can make in answering this question, which is often the case when translating a problem statement into a simulation. In our first tutorial, we will approximate the problem by the following:

• The (initial) position and velocity of Asterix are given in the form of a Keplerian state vector.
• The Earth is modeled as a point mass, and the spacecraft is modeled as massless.
• All perturbations, such as solar radiation and third-body gravity, etc. are neglected.

Note

In a followup tutorial, we will show how to extend the dynamical model of the spacecraft, slightly relaxing the last assumption.

The main focus of this tutorial is to provide insight into the use of the Tudat library for your application. Using this test case, we will show how you can use some of the elements of the Tudat library, namely:

• The transformation from Keplerian elements (provided) to Cartesian elements (required by propagator).
• The settings for the environment
• The settings for the acceleration models.
• The settings for the state derivative model (propagator).
• The settings for the numerical integrator.
• The use of ‘create’ functions to parse these settings

In particular, it will show you how these elements can work together to provide you with a powerful satellite propagator tool. It is not the intention to give you a C++ tutorial. Only the parts of the source code specific to Tudat will be discussed on this page. Detailed documentation on the environment models, acceleration models and propagators/integrators is available, but we recommend working through the tutorials here before delving into the details of all options that are available to you.

In Tudat, the environment, acceleration models and propagator are not created directly by the user. Instead, the user defines a number of ‘Settings’ objects that describe how the actual models are to be put together. These Settings objects are then passed to ‘create’ functions to build up the objects/data used in the simulations, and link them all together. In linking these objects together, information is automatically updated, and where necessary transformed and processed in the manner specified by the user in the Settings objects.

1.1. Set up the environment

First we discuss the setup of the environment, which is stored in a list of Body objects, each of which represents either a natural celestial body, or a spacecraft (orbiter, entry vehicle, etc.). The entire environment (gravity fields, atmospheres, ephemerides, rotation models, etc.) are then defined as members of the corresponding Body object. For this example, the only environment models we need are:

• Earth gravity field (point-mass)
• Earth ephemeris model

The block of code used to generate the environment that we require is:

// Create body objects.
std::vector< std::string > bodiesToCreate;
bodiesToCreate.push_back( "Earth" );
std::map< std::string, boost::shared_ptr< BodySettings > > bodySettings =
        getDefaultBodySettings( bodiesToCreate );
bodySettings[ "Earth" ]->ephemerisSettings = boost::make_shared< ConstantEphemerisSettings >(
            Eigen::Vector6d::Zero( ) );

// Create Earth object
NamedBodyMap bodyMap = createBodies( bodySettings );

// Create spacecraft object.
bodyMap[ "Asterix" ] = boost::make_shared< simulation_setup::Body >( );

// Finalize body creation.
setGlobalFrameBodyEphemerides( bodyMap, "SSB", "J2000" );

Creating an environment starts by creating the settings for all required environment models, which are stored in a list of BodySettings objects. Each BodySettings object has a list of objects for environment model settings, which are empty upon creation of the BodySettings, and may be set by the user to their desired specifications. To save you the trouble of having to define all settings manually, we provide default options for the environment, many of which will often be sufficient for a given simulation. A list of all possible properties of the BodySettings, as well as the default settings, can be found here. In our example, we only need environment models for the Earth, so we start by:

std::vector< std::string > bodiesToCreate;
bodiesToCreate.push_back( "Earth" );
std::map< std::string, boost::shared_ptr< BodySettings > > bodySettings =
        getDefaultBodySettings( bodiesToCreate );

With the settings that are now stored in the bodySettings map, we could generate our environment and move on to the next step of the simulation. However, both to showcase one of the options of setting up the environment, and for numerical accuracy, we make one modification to the default ephemeris settings, defining the Earth to be fixed at the center of the Solar system:

bodySettings[ "Earth" ]->ephemerisSettings = boost::make_shared< ConstantEphemerisSettings >(
        Eigen::Vector6d::Zero( ), "SSB", "J2000" );

Warning

By using this code, we ‘cheat’ a little bit, since we put the Earth in the main inertial frame for Tudat: the Solar System Barycenter. This approach should only be used when considering no third-body perturbations. Note that we could have used any ephemeris setting for the Earth (since we propagate our satellite w.r.t. the Earth origin) without changing anything in our dynamical model, but at a slight loss of numerical precision.

Now, we have created the settings we need for the environment, and we can move on to creating the environment models themselves, which are stored in a set of Body objects. These Body objects are stored in a NamedBodyMap, which is a typedef (shorthand name) for:

std::unordered_map< std::string, boost::shared_ptr< simulation_setup::Body > >

This unordered_map may be accessed as a regular map. The std::string keys represent the names of the bodies in the list, and the value boost::shared_ptr the corresponding Body object containing all the environment models. With the settings of our ephemeris and gravity field, we now create the bodyMap by:

NamedBodyMap bodyMap = createBodies( bodySettings );

1.2. Create the vehicle

Our environment is now missing only one aspect: the spacecraft. Our spacecraft (called Asterix) requires no specific properties, it merely needs to exist (its initial state is defined in a subsequent part of the example). Therefore, we can simply add the Asterix satellite to our environment by creating a new empty Body object:

bodyMap[ "Asterix" ] = boost::make_shared< simulation_setup::Body >( );

Although not required in this simulation, it is good practice to call the following function following the complete setup of the bodyMap:

setGlobalFrameBodyEphemerides( bodyMap, "SSB", "J2000" );

Calling this function will allow hierarchical ephemerides to be properly used in the simulation (i.e. orbiter ephemeris w.r.t. Moon, Moon w.r.t. Earth, Earth w.r.t. Sun, Sun w.r.t. barycenter).

1.3. Set up the acceleration models

To define the settings of the propagation of the orbit, we start by defining the required acceleration models. The block of code that performs the required operations is:

// Define propagator settings variables.
SelectedAccelerationMap accelerationMap;
std::vector< std::string > bodiesToPropagate;
std::vector< std::string > centralBodies;

bodiesToPropagate.push_back( "Asterix" );
centralBodies.push_back( "Earth" );

// Define propagation settings.
std::map< std::string, std::vector< boost::shared_ptr< AccelerationSettings > > > accelerationsOfAsterix;
accelerationsOfAsterix[ "Earth" ].push_back( boost::make_shared< AccelerationSettings >(
                                                 basic_astrodynamics::central_gravity ) );
accelerationMap[  "Asterix" ] = accelerationsOfAsterix;

// Create acceleration models and propagation settings.
basic_astrodynamics::AccelerationMap accelerationModelMap = createAccelerationModelsMap(
            bodyMap, accelerationMap, bodiesToPropagate, centralBodies );

As with the environment models, there is no need to manually create the models. The user must only define the properties of the acceleration models that are desired, which are:

• List of bodies that are to be numerically propagated.
• Origin of reference frame in which they are to be propagated (may be different for each body). Note that propagation in Tudat is always done in a non-rotating reference frame, only the origin of the frames can be varied.
• A list of settings of the accelerations model(s) acting on the bodies in the simulation.

These properties are to be defined in the following variables:

// Define propagator settings variables.
SelectedAccelerationMap accelerationMap;
std::vector< std::string > bodiesToPropagate;
std::vector< std::string > centralBodies;

The list of propagated bodies and the reference frame origins (central bodies) are simply lists of strings, which for our case are:

bodiesToPropagate.push_back( "Asterix" );
centralBodies.push_back( "Earth" );

These settings for the accelerations require some more structure, though, and are stored in a SelectedAccelerationMap. This is typedef for:

std::map< std::string, std::map< std::string, std::vector< boost::shared_ptr< AccelerationSettings > > > >

This is a double map (with twice a string as a key). The two levels correspond to the names of bodies undergoing an acceleration (first key), and those for bodies exerting an acceleration (second key). This allows any number of bodies to be propagated, undergoing any number (and type) of accelerations. Mutual acceleration between bodies being propagated, as is the case for Solar system dynamics for instance, is automatically handled by the code and requires no specific consideration.

In our example, we have only a single point-mass acceleration due to Earth, acting on Asterix. We define the settings for the acceleration as follows:

std::map< std::string, std::vector< boost::shared_ptr< AccelerationSettings > > > accelerationsOfAsterix;
accelerationsOfAsterix[ "Earth" ].push_back( boost::make_shared< AccelerationSettings >(
                                                 basic_astrodynamics::central_gravity ) );
accelerationMap[  "Asterix" ] = accelerationsOfAsterix;

A single acceleration, of type central_gravity to be exerted by body "Earth" on body "Asterix" is now defined. The list of the actual acceleration models is now created by:

basic_astrodynamics::AccelerationMap accelerationModelMap = createAccelerationModelsMap(
            bodyMap, accelerationMap, bodiesToPropagate, centralBodies );

which automatically links together all required objects and functions.

1.4. Set up the propagation settings

Now that we have both our environment models and our acceleration model, we can create the full settings for the propagation. These settings are stored in a PropagatorSettings object. For this example, we will only consider the propagation of translational dynamics, which is stored in the derived class TranslationalStatePropagatorSettings. The settings for the propagator are the following:

• The acceleration models.
• The list of bodies that are to be propagated.
• The origins w.r.t. which these bodies are to be propagated.
• The initial Cartesian state that is to be used.
• Termination conditions for the propagation (here, a fixed final time).

The above settings are provided in the following block of code:

// Set Keplerian elements for Asterix.
Vector6d asterixInitialStateInKeplerianElements;
asterixInitialStateInKeplerianElements( semiMajorAxisIndex ) = 7500.0E3;
asterixInitialStateInKeplerianElements( eccentricityIndex ) = 0.1;
asterixInitialStateInKeplerianElements( inclinationIndex ) = convertDegreesToRadians( 85.3 );
asterixInitialStateInKeplerianElements( argumentOfPeriapsisIndex )
        = convertDegreesToRadians( 235.7 );
asterixInitialStateInKeplerianElements( longitudeOfAscendingNodeIndex )
        = convertDegreesToRadians( 23.4 );
asterixInitialStateInKeplerianElements( trueAnomalyIndex ) = convertDegreesToRadians( 139.87 );

// Convert Asterix state from Keplerian elements to Cartesian elements.
double earthGravitationalParameter = bodyMap.at( "Earth" )->getGravityFieldModel( )->getGravitationalParameter( );
Eigen::VectorXd systemInitialState = convertKeplerianToCartesianElements(
            asterixInitialStateInKeplerianElements,
            earthGravitationalParameter );

boost::shared_ptr< TranslationalStatePropagatorSettings< double > > propagatorSettings =
        boost::make_shared< TranslationalStatePropagatorSettings< double > >
        ( centralBodies, accelerationModelMap, bodiesToPropagate, systemInitialState, simulationEndEpoch );

If the body that is being propagated has a pre-existing ephemeris, the initial state may be retrieved automatically. In this example, however, we manually define our initial state from the Keplerian state:

// Set Keplerian elements for Asterix.
Vector6d asterixInitialStateInKeplerianElements;
asterixInitialStateInKeplerianElements( semiMajorAxisIndex ) = 7500.0E3;
asterixInitialStateInKeplerianElements( eccentricityIndex ) = 0.1;
asterixInitialStateInKeplerianElements( inclinationIndex ) = convertDegreesToRadians( 85.3 );
asterixInitialStateInKeplerianElements( argumentOfPeriapsisIndex )
        = convertDegreesToRadians( 235.7 );
asterixInitialStateInKeplerianElements( longitudeOfAscendingNodeIndex )
        = convertDegreesToRadians( 23.4 );
asterixInitialStateInKeplerianElements( trueAnomalyIndex ) = convertDegreesToRadians( 139.87 );

// Convert Asterix state from Keplerian elements to Cartesian elements.
double earthGravitationalParameter = bodyMap.at( "Earth" )->getGravityFieldModel( )->getGravitationalParameter( );
Eigen::VectorXd systemInitialState = convertKeplerianToCartesianElements(
            asterixInitialStateInKeplerianElements,
            earthGravitationalParameter );

Note that we use the Earth gravity field inside the bodyMap for the conversion from Keplerian to Cartesian coordinates. Now, we can create our propagator settings by:

boost::shared_ptr< TranslationalStatePropagatorSettings< > > propagatorSettings =
    boost::make_shared< TranslationalStatePropagatorSettings< > >
        ( centralBodies, accelerationModelMap, bodiesToPropagate, systemInitialState, simulationEndEpoch);

Where we have passed exactly the five aspects listed above as input to the TranslationalStatePropagatorSettings. If you have a look at the code for the TranslationalStatePropagatorSettings, you will notice that there are multiple constructors for the class, each with a number of additional input arguments (for which we use the default values). These more advanced options are discussed in the following tutorials.

A final piece of information needed to propagate the orbit is the settings object for the numerical integration. We use a Runge-Kutta 4 integrator, with a 10 second time step, starting the numerical integration at t = 0:

// Create numerical integrator.
double simulationStartEpoch = 0.0;
const double fixedStepSize = 10.0;
boost::shared_ptr< IntegratorSettings< > > integratorSettings =
    boost::make_shared< IntegratorSettings< > >
        ( rungeKutta4, simulationStartEpoch, fixedStepSize );

1.5. Perform the orbit propagation

Now, we have defined all the information needed to propagate the orbit of our satellite, which are stored in the bodyMap (environment), propagatorSettings (settings for the full state derivative model) and integratorSettings (settings on how to obtain the numerical solution). The propagation is done by an object of type (derived from) DynamicsSimulator:

SingleArcDynamicsSimulator< > dynamicsSimulator(
        bodyMap, integratorSettings, propagatorSettings );

Upon creating this class, the numerical propagation is performed, and the output is stored in the class. Various options exist for parsing the output of the numerical propagation, which will be discussed in the next tutorials. The numerical solution of the orbit can be retrieved as follows:

std::map< double, Eigen::VectorXd > integrationResult = dynamicsSimulator.getEquationsOfMotionNumericalSolution( );

where the retrieved result is a std::map where the key double is the time at each step in the integration, and the value, Eigen::VectorXd, is the corresponding Cartesian state of Asterix w.r.t. the Earth, in the J2000 reference frame. To analyze/plot your numerical results further using e.g. Matlab, you can print the output to a text file as follows:

// Write satellite propagation history to file.
input_output::writeDataMapToTextFile( integrationResult,
                                      "singleSatellitePropagationHistory.dat",
                                      tudat_applications::getOutputPath( ),
                                      "",
                                      std::numeric_limits< double >::digits10,
                                      std::numeric_limits< double >::digits10,
                                      "," );

1.6. Results

The output of the program should look similar to the output below:

Starting ../tudatBundle/tudatExampleApplications/satellitePropagatorExamples/bin/applications/application_SingleSatellitePropagator...
Single Earth-Orbiting Satellite Example.
The initial position vector of Asterix is [km]:
7037.48
3238.06
2150.72
The initial velocity vector of Asterix is [km/s]:
  -1.46566
-0.0409584
    6.6228
After 86400 seconds, the position vector of Asterix is [km]:
-4560.45
-1438.32
 5973.99
And the velocity vector of Asterix is [km/s]:
-4.55021
-2.41254
-4.95063
../tudatBundle/tudatExampleApplications/satellitePropagatorExamples/bin/applications/application_SingleSatellitePropagator exited with code 0