4. Inner Solar System Propagation

The example described on this page aims to simulate the dynamics of the main celestial bodies in the inner solar system. The code for this tutorial is located in your Tudat Bundle at:


For this example, we have the following problem statement:

Given the position and velocity of the main celestial bodies in the inner solar system at a certain point in time, what will their position and velocity after a certain period of time?


The example described in this page assumes that the user has read the Unperturbed Earth-orbiting Satellite. This page only describes the differences with respect to such example, so please go back before proceding.

4.1. Set up the environment

The first step is to define the bodies that will be propagated, which are stored hierarchically in a moons to stars fashion. The reason why this is done will become clear later:

// Define bodies in simulation.
unsigned int totalNumberOfBodies = 6;
std::vector< std::string > bodyNames;
bodyNames.resize( totalNumberOfBodies );
bodyNames[ 0 ] = "Moon";
bodyNames[ 1 ] = "Earth";
bodyNames[ 2 ] = "Mars";
bodyNames[ 3 ] = "Venus";
bodyNames[ 4 ] = "Mercury";
bodyNames[ 5 ] = "Sun";

The next step is to save such bodies in the bodyMap and set the global ephemerides:

// Create bodies needed in simulation
std::map< std::string, boost::shared_ptr< BodySettings > > bodySettings =
        getDefaultBodySettings( bodyNames );
NamedBodyMap bodyMap = createBodies( bodySettings );
setGlobalFrameBodyEphemerides( bodyMap, "SSB", "ECLIPJ2000" );

4.2. Set up the acceleration models

The example explained in this page is slightly different than all the previous examples, since two simulations are performed within the same application:

  1. The Solar System Barycenter (SSB) is considered as the central body for all the propagated celestial bodies.
  2. The propagated celestial bodies consider the higher hierarchical body as their respective central body. (Moon > Earth; Mercury, Venus, Earth, Mars > Sun; Sun > SSB).

These two simulations result in two different output files. Consequently, all the steps required to create a DynamicsSimulator are performed twice, as shown by the following for loop:

// Run simulation for 2 different central body settings (barycentric and hierarchical)
for( int centralBodySettings = 0; centralBodySettings < 2; centralBodySettings++ )
    ///////////////////////             CREATE ACCELERATIONS            ///////////////////////////////////////////////

    // Set accelerations between bodies that are to be taken into account (mutual point mass gravity between all bodies).
    SelectedAccelerationMap accelerationMap;
    for( unsigned int i = 0; i < bodyNames.size( ); i++ )
        std::map< std::string, std::vector< boost::shared_ptr< AccelerationSettings > > > currentAccelerations;
        for( unsigned int j = 0; j < bodyNames.size( ); j++ )
            // Create central gravity acceleration between each 2 bodies.
            if( i != j )
                currentAccelerations[ bodyNames.at( j ) ].push_back(
                            boost::make_shared< AccelerationSettings >( central_gravity ) );\
        accelerationMap[ bodyNames.at( i ) ] = currentAccelerations;

The first step in the for loop is to create the acceleration models. The only acceleration model used in this example is the central_gravity model, which in this case acts on all bodies from all bodies. The if condition shown above ensures that the gravity field model of a body does not act on itself.


The reason behind having two different simulations is merely for illustration purposes and is not requirement when simulating the Solar System dynamics.

4.3. Set up the propagation settings

The next step found in the for loop is to create the list with the bodiesToPropagate. Contrary to the Perturbed Earth-orbiting Satellite example, in this case all the bodies in the bodyMap are propagated:

// Define list of bodies to propagate
std::vector< std::string > bodiesToPropagate = bodyNames;
unsigned int numberOfNumericalBodies = bodiesToPropagate.size( );

Since all bodies require a central body, a placeholder is created for each body:

// Define central bodies to use in propagation.
std::vector< std::string > centralBodies;
centralBodies.resize( numberOfNumericalBodies );

where each place holder is filled depending on which simulation is being executed. The first part of the following if loop corresponds to Simulation 1 and the second part (the elseif loop) corresponds to Simulation 2:

// Set central body as Solar System Barycente for each body
if( centralBodySettings == 0 )
   for( unsigned int i = 0; i < numberOfNumericalBodies; i++ )
      centralBodies[ i ] = "SSB";
else if( centralBodySettings == 1 )
   for( unsigned int i = 0; i < numberOfNumericalBodies; i++ )
      // Set Earth as central body for Moon
      if( i == 0 )
         centralBodies[ i ] = "Earth";
      // Set barycenter as central 'body' for Sun
      else if( i == 5 )
         centralBodies[ i ] = "SSB";
      // Set Sun as central body for all planets
         centralBodies[ i ] = "Sun";

Once the centralBodies have been set, the creation of the AccelerationMap can be finalized:

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

The next step in the for loop is to define the propagation time. In this example, the orbits of the inner solar system bodies are initiated 1.0e7 seconds after J2000 and are propagated for five Julian years:

// Specify initial time
double initialEphemerisTime = 1.0E7;
double finalEphemerisTime = 1.0E7 + 5.0 * physical_constants::JULIAN_YEAR;

// Get initial state vector as input to integration.
Eigen::VectorXd systemInitialState = getInitialStatesOfBodies(
   bodiesToPropagate, centralBodies, bodyMap, initialEphemerisTime );

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

Once this step is completed, the DynamicsSimulator object is created and the equations of motion are integrated. Once again, this step is within the for loop, thus this step is performed twice:

// Create simulation object and propagate dynamics.
SingleArcDynamicsSimulator< > dynamicsSimulator(
   bodyMap, integratorSettings, propagatorSettings, true, false, false );

4.3.1. Save the propagation results

The final step in the for loop is to save the propagation results to a file. First, the propagation history is retrieved from the dynamicsSimulator:

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

where the retrieved result is an std::map with a double as key (the simulation time) and a Eigen::VectorXd of 1x36 elements, containing the 6-variable state for the six bodies. Next, an std::vector where each of its elements is an std::map that contains the propagation history of a single body:

// Retrieve numerically integrated state for each body.
std::vector< std::map< double, Eigen::VectorXd > > allBodiesPropagationHistory;
allBodiesPropagationHistory.resize( numberOfNumericalBodies );

for( std::map< double, Eigen::VectorXd >::const_iterator stateIterator = integrationResult.begin( );
      stateIterator != integrationResult.end( ); stateIterator++ )
   for( unsigned int i = 0; i < numberOfNumericalBodies; i++ )
      allBodiesPropagationHistory[ i ][ stateIterator->first ] = stateIterator->second.segment( i * 6, 6 );

Finally, each of the std::map within allBodiesPropagationHistory are written to individual .dat files:

 for( unsigned int i = 0; i < numberOfNumericalBodies; i++ )
    // Write propagation history to file.
       allBodiesPropagationHistory[ i ],
       "innerSolarSystemPropagationHistory" + bodyNames.at( i ) +
       boost::lexical_cast< std::string >( centralBodySettings ) + ".dat",
       tudat_applications::getOutputPath( ),
       std::numeric_limits< double >::digits10,
       std::numeric_limits< double >::digits10,
       "," );