In some cases, the various features of Tudat do not contain the tool(s) to solve a specific problem. As Tudat is written in C++ and all the source code is readily available to the user, it is possible to add features to the Tudat software. This does require a good understanding of the C++ language and of how Tudat works, thus it is only recommended to users who have been working with Tudat before.

## Adding Dependent Variables to Save¶

As the execution of the simulator happens “behind closed doors” in the DynamicsSimulator function, it can be non-trivial to retrieve the history of dependent variables. Luckily, as explained in Propagator Settings: Dependent Variables, it is possible to create a list of dependent variables and use it as an input to the propagation settings to be able to retrieve a history of these dependent variables after the simulation is done. A list of the available dependent variables is given in: Propagator Settings: Dependent Variables. It is possible that another dependent variable is needed by the user, which is not given in this list. This guide will show how a new dependent variable can be added to the code.

To add a new dependent variable, four files will need to be changed, these are all located in the following directory:

…/tudatBundle/tudat/Tudat/SimulationSetup/PropagationSetup/

For this guide, the mach_number_dependent_variable will be used as an example of how a dependent variable is implemented. Be aware that there are some slight differences between the implementation of the mach_number_dependent_variable and other variables, but this will be explained in the guide.

The first addition needs to be made in the propagationOutputSettings.h file. Here, an enum called PropagationDependentVariables is given which shows the names of all the available dependent variables that can be saved. The list looks as follows;

enum PropagationDependentVariables
{
mach_number_dependent_variable = 0,
altitude_dependent_variable = 1,
airspeed_dependent_variable = 2,
...
...
...
body_fixed_groundspeed_based_velocity_variable = 31,
keplerian_state_dependent_variable = 32,
modified_equinocial_state_dependent_variable = 33
};


The name of the new variable can thus be added to the end of the list, with a number assigned to it.

The next addition needs to be made in the propagationOutputSettings.cpp file. In this file, a method is present called: getDependentVariableName. This method contains a switch statement which contains cases for all the different dependent variables. For a new variable, a new case must be made with the name that was added to the PropagationDependentVariables list. Inside the case, a string called variableName must be assigned with the name of the to be added variable. In the case of the mach number, it looks as follows:

case mach_number_dependent_variable:
variableName = "Mach number ";
break;


It is important to remember that the case should be the same as the name given in PropagationDependentVariables, but the variableName can be chosen to the preference of the user.

The third addition needs to be made in propagationOutput.cpp. Here a method called getDependentVariableSize is located. In this function, another switch statement is made with the same cases as before. Again a new case should be made for the new variable, but now, inside the case, another variable called variableSize should be changed to the size of the added variable. Thus this could be 1 for a double, or 3 for an Eigen::Vector3d. In the case of the mach number, it looks as follows:

case mach_number_dependent_variable:
variableSize = 1;
break;


The final change should then be made to the propagationOutput.h file. In a method called getDoubleDependentVariableFunction, the actual calculation of the variable is done. Again, a switch stamentent in the same manner as before is made, with in every case the calculation of the variable. This method gets the bodyMap as input, thus methods available inside the FlightConditions could, for example, be used to calculate the dependent variable. For the new variable, a new case needs to be made in which a std::function< double( )> is returned. This function can take several variables as input and should return the dependent variable. The implementation for the mach number is given here:

  case mach_number_dependent_variable:
{
if( bodyMap.at( bodyWithProperty )->getFlightConditions( ) == NULL )
{
std::string errorMessage = "Error, no flight conditions available when requesting Mach number output of " +
bodyWithProperty + "w.r.t." + secondaryBody;
throw std::runtime_error( errorMessage );
}

std::function< double( const double, const double ) > functionToEvaluate =
std::bind( &aerodynamics::computeMachNumber, std::placeholders::_1, std::placeholders::_2 );

// Retrieve functions for airspeed and speed of sound.
std::function< double( ) > firstInput =
std::bind( &aerodynamics::FlightConditions::getCurrentAirspeed,
bodyMap.at( bodyWithProperty )->getFlightConditions( ) );
std::function< double( ) > secondInput =
std::bind( &aerodynamics::FlightConditions::getCurrentSpeedOfSound,
bodyMap.at( bodyWithProperty )->getFlightConditions( ) );

variableFunction = std::bind( &evaluateBivariateFunction< double, double >,
functionToEvaluate, firstInput, secondInput );
break;
}


If the variable is a vector (or matrix) and not a double, the case should not be added to getDoubleDependentVariableFunction, but to: getVectorDependentVariableFunction

If this is all done, the dependent variable name can be added to the dependent variables save list.

## Adding a New Environment Model¶

Tudat contains various pre-built environment models that users are able to implement in their simulations. However, for some applications, it might be necesarry to alter an exisiting environment model, or add a completely new one. This section will describe how a new environment model can be added to Tudat.

First, it is important to understand how environment models are used in Tudat, this will give some insight in what needs to be added/changed if a new environment model is integrated into Tudat. For this tutorial, an existing environment model will be used as an example: the tabulated atmosphere model.

A user will define a list of bodies to be used in the simulation, an example is given here:

// Define simulation body settings.
std::map< std::string, std::shared_ptr< BodySettings > > bodySettings =
getDefaultBodySettings( { "Mars" }, simulationStartEpoch - 10.0 * fixedStepSize,
simulationEndEpoch + 10.0 * fixedStepSize );
bodySettings[ "Mars" ]->ephemerisSettings = std::make_shared< simulation_setup::ConstantEphemerisSettings >(
Eigen::Vector6d::Zero( ), "SSB", "J2000" );
bodySettings[ "Mars" ]->rotationModelSettings->resetOriginalFrame( "J2000" );
std::string atmosphereFile = ... ;
bodySettings[ "Mars" ]->atmosphereSettings = std::make_shared< TabulatedAtmosphereSettings >( atmosphereFile );


The BodySettings class contains a member variable called atmosphereSettings of type AtmosphereSettings. As can be seen in the last line of the example code above, this member variable is set to a TabulatedAtmosphereSettings by the user. This is a valid statement (eventhough atmosphereSettings is not of type TabulatedAtmosphereSettings) because TabulatedAtmosphereSettings is derived from the base class AtmosphereSettings. These concepts are called polymorphism and inheritance, and should be understood by the reader before continuing.

When all the bodySettings are defined by the user, the following command will be executed:

simulation_setup::NamedBodyMap bodyMap = simulation_setup::createBodies( bodySettings );


The simulation_setup::createBodies() will set all the environment models up to be used later in the simulation, and stores them in the variable bodyMap. If the simulation_setup::createBodies() is looked at, the following can be seen:

// Create atmosphere model objects for each body (if required).
for( unsigned int i = 0; i < orderedBodySettings.size( ); i++ )
{
if( orderedBodySettings.at( i ).second->atmosphereSettings != NULL )
{
bodyMap[ orderedBodySettings.at( i ).first ]->setAtmosphereModel(
createAtmosphereModel( orderedBodySettings.at( i ).second->atmosphereSettings,
orderedBodySettings.at( i ).first ) );
}
}


Where orderedBodySettings is an ordered map of all the bodies. This code goes through all the bodies and checks if the bodies have an atmosphere model set or not. If not, the atmosphere model is set by the setAtmosphereModel() function, which takes as input a pointer to the atmosphereModel of the body. This pointer is returned by the function: createAtmosphereModel(), which creates the atmosphere model using the settings defined by the user in the variable atmosphereSettings. The createAtmosphereModel() looks as follows:

std::shared_ptr< aerodynamics::AtmosphereModel > createAtmosphereModel(
const std::shared_ptr< AtmosphereSettings > atmosphereSettings,
const std::string& body )
{
using namespace tudat::aerodynamics;

// Declare return object.
std::shared_ptr< AtmosphereModel > atmosphereModel;

// Check which type of atmosphere model is to be created.
switch( atmosphereSettings->getAtmosphereType( ) )
{
case exponential_atmosphere:
{
...
...
...
}
case tabulated_atmosphere:
{
// Check whether settings for atmosphere are consistent with its type
std::shared_ptr< TabulatedAtmosphereSettings > tabulatedAtmosphereSettings =
std::dynamic_pointer_cast< TabulatedAtmosphereSettings >( atmosphereSettings );
if( tabulatedAtmosphereSettings == NULL )
{
throw std::runtime_error(
"Error, expected tabulated atmosphere settings for body " + body );
}
else
{
// Create and initialize tabulatedl atmosphere model.
atmosphereModel = std::make_shared< TabulatedAtmosphere >(
tabulatedAtmosphereSettings->getAtmosphereFile( ) );
}
break;
}
case nrlmsise00:
{
...
...
...
}
default:
throw std::runtime_error(
"Error, did not recognize atmosphere model settings type " +
std::to_string( atmosphereSettings->getAtmosphereType( ) ) );
}

return atmosphereModel;
}


This function checks which atmosphere model is used by using a switch statement with the atmosphereSettings->getAtmosphereType( ). The three cases are defined using the enum AtmosphereTypes (located in createAtmosphereModel.h). In the case of the tabulated atmosphere, first a check is made if tabulated atmosphere setting are actually initialized. If they are, the atmosphereModel is set to a tabulated atmosphere using the settings from atmosphereSettings. It can be seen in the example code above that the variable atmosphereModel is of type AtmosphereModel, but it is set to a variable of type TabulatedAtmosphere in the tabulated atmosphere case. This is again valid due to the fact that the class TabulatedAtmosphere is derived from the base class AtmosphereModel.

After the atmosphere model is set in the body map, it can be used whenever a certain quantity, e.g. the density or temperature, is needed by another part of the simulation. How these quantities are calculated is defined in the TabulatedAtmosphere class.

Now that the creation of an environment model is understood, it can be discussed what should change in the Tudat code when something is added. There are two options when changing the code: the user can modify an existing environment model, or they can add a new model to an existing type of environment model. These type of modifications require different amount of changes made to Tudat and are thus explained seperately here:

• Modify existing environment model: this option is the easiest as it only requires changes in the specific environment type files. If a function is added to an environment model, it is important to also include this function in the base class file of that specific environment model, with the virtual statement. For example, take the getDensity() function. This function is put in the file AtmosphereModel.h as a pure virtual function (don’t need a function definition), by putting a =0 after the function definition and including the virtual statement. Now, in every derived class, this function should return something, depending on the implementation. If only something inside an already existing function needs to be changed, it (most of the time) shouldn’t be changed in all the derived classes. If a variable is added to the constructor of the class, it is important that all the cases that this class is called should be changed in the entire Tudat code to prevent errors (use the find usages option in Qt to find them).
• Add a new environment model: this modification requires some extra changes to the framework of the environment model implementation. First, when the new model is made, make sure that it is derived from the base class, and that it contains some of the basic functions. Once the new model is made, a Settings class should be made in the same way as TabulatedAtmosphereSettings. This class should be added to the create... files, and should contain functions that store variables that are used to call the constructor of the corresponding class (again, make sure it is derived from the proper base class). Then, in the corresponding .cpp file, a new case for the new environment model should be added to the create... function, just as in the tabulated atmosphere example. Make sure that in the getAtmosphereType( ) function, the name of the new model is also included. In the case statement, make sure to add checks, and throw runtime errors if they are violated. Another step that needs to be taken is to update the createEnvironmentUpdater.cpp file. This file includes several switch statements that need to have the new acceleration model in it. Use the existing code to determine how the new case should be made.

## Adding a New Acceleration Model¶

For some applications, the various acceleration models in Tudat might not be sufficient. This guide will explain how acceleration models are set-up within Tudat, and how these models can be modified, or how to add a completely new model.

Just as before, it is important to understand the framework of the acceleration models. This will be explained using an example, namely the aerodynamic acceleration model. The first time the acceleration models will be needed is in the user’s own code. An example of how this can be done is shown below:

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

// Define acceleration model settings.
std::map< std::string, std::vector< std::shared_ptr< AccelerationSettings > > > accelerationsOfWaverider;
accelerationsOfWaverider[ "Mars" ].push_back( std::make_shared< AccelerationSettings >( aerodynamic ) );
accelerationMap[  "Vehicle" ] = accelerationsOfWaverider;

bodiesToPropagate.push_back( "Vehicle" );
centralBodies.push_back( "Mars" );

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


The user will define a special map, with type SelectedAccelerationMap that works as follows: the first key defines the body exerting the acceleration, the second key the body undergoing the acceleration, and the value in the map is a vector of pointers to a type class called AccelerationSettings, which takes as input an enum called AvailableAcceleration, which lists all the available acceleration models. This map is then used as an input for the createAccelerationModelsMap function, which is a function that creates a special acceleration map that can be used as an input to the propagator settings.

There are two functions, located in CreateAccelerationModels.cpp, called createAccelerationModelsMap. The function that is called will first order the map, and then return the value of the other function, which takes as input the ordered map, made in the first createAccelerationModelsMap. The second function looks as follows:

//! Function to create a set of acceleration models from a map of bodies and acceleration model types.
basic_astrodynamics::AccelerationMap createAccelerationModelsMap(
const NamedBodyMap& bodyMap,
const SelectedAccelerationMap& selectedAccelerationPerBody,
const std::map< std::string, std::string >& centralBodies )
{
...
...
...
currentAcceleration = createAccelerationModel( bodyMap.at( bodyUndergoingAcceleration ),
bodyMap.at( bodyExertingAcceleration ),
accelerationsForBody.at( i ).second,
bodyUndergoingAcceleration,
bodyExertingAcceleration,
currentCentralBody,
currentCentralBodyName,
bodyMap );

// Create acceleration model.
mapOfAccelerationsForBody[ bodyExertingAcceleration ].push_back(
currentAcceleration );

...
...
...
// Put acceleration models on current body in return map.
accelerationModelMap[ bodyUndergoingAcceleration ] = mapOfAccelerationsForBody;
}
return accelerationModelMap;
}


This function does several checks if the right models are present, and then calls the createAccelerationModel( ) function. This function is where the acceleration models are called and looks as follows:

//! Function to create acceleration model object.
std::shared_ptr< AccelerationModel< Eigen::Vector3d > > createAccelerationModel(
const std::shared_ptr< Body > bodyUndergoingAcceleration,
const std::shared_ptr< Body > bodyExertingAcceleration,
const std::shared_ptr< AccelerationSettings > accelerationSettings,
const std::string& nameOfBodyUndergoingAcceleration,
const std::string& nameOfBodyExertingAcceleration,
const std::shared_ptr< Body > centralBody,
const std::string& nameOfCentralBody,
const NamedBodyMap& bodyMap )
{
// Declare pointer to return object.
std::shared_ptr< AccelerationModel< Eigen::Vector3d > > accelerationModelPointer;

// Switch to call correct acceleration model type factory function.
switch( accelerationSettings->accelerationType_ )
{
...
...
...
case aerodynamic:
accelerationModelPointer = createAerodynamicAcceleratioModel(
bodyUndergoingAcceleration,
bodyExertingAcceleration,
nameOfBodyUndergoingAcceleration,
nameOfBodyExertingAcceleration );
break;
...
...
...
default:
throw std::runtime_error(
std::string( "Error, acceleration model ") +
std::to_string( accelerationSettings->accelerationType_ ) +
" not recognized when making acceleration model of" +
nameOfBodyExertingAcceleration + " on " +
nameOfBodyUndergoingAcceleration );
break;
}
return accelerationModelPointer;
}


Ths function uses a switch statement on the accelerationSettings->accelerationType_ (which was set by the user in the beginning), to determine which acceleration model needs to be called. The createAerodynamicAcceleratioModel( ) then does several checks to determine if all the models that are used to calculate the aerodynamic force are present, after which it calls the constructor of the AerodynamicAcceleration class using the input values created in createAerodynamicAcceleratioModel( ). The structure of the AerodynamicAcceleration class, and how one can make a new acceleration model will be discussed in the next section.

When making a new acceleration model, there are two things that need to be added first. The first one is a new class containing functions to calculate the accelerations acting on the respective body. It is important that this class is derived from the base class: basic_astrodynamics::AccelerationModel< dataType >, where dataType is the type of the acceleration variable. Furthermore, there are two functions that need to be present in the new acceleration model class:

• void updateMembers(const double currentTime = TUDAT_NAN). This function updates all the member variables to the current situation so they can be used by the other necessary function.
• dataType getAcceleration( ). This function uses the recently updated member variables to calculate the acceleration acting on the body.

For the case of the aerodynamic acceleration, these two functions look as follows:

Eigen::Vector3d getAcceleration( )
{
return computeAerodynamicAcceleration(
0.5 * currentDensity_ * currentAirspeed_ * currentAirspeed_,
currentReferenceArea_, currentForceCoefficients_, currentMass_ );
}

void updateMembers( const double currentTime = TUDAT_NAN )
{
if( !( this->currentTime_ == currentTime ) )
{
currentForceCoefficients_ = coefficientMultiplier_ * this->coefficientFunction_( );
currentDensity_ = this->densityFunction_( );
currentMass_ = this->massFunction_( );
currentAirspeed_ = this->airSpeedFunction_( );
currentReferenceArea_ = this->referenceAreaFunction_( );

currentTime_ = currentTime;
}
}


where computeAcceleration() is a function that uses the well known lift, side, and drag acceleration equations to calculate the acceleration vector.

The second piece of code that needs to be added is the new create function inside createAccelerationModel.cpp. This function should return a pointer to the respective acceleration class, and takes as input the various models needed for the calculation of the acceleration. The main goal of this function is to do some checks if the necessary models are present, and either throw an error, or initialize the model there. When this is done, it uses all the available models as input to the acceleration model, and returns a pointer to this class.

Once the acceleration class is made, and the create function is in place, the acceleration model should be implemented into the acceleration framework. The first step for this is to add the a new model to the AvailableAcceleration enum. This name could be anything, as long as it relates to the actual acceleration model. Nothing has to be added to the createAccelerationModelsMap, however, the createAccelerationModel function has a switch statement that needs to be altered. This switch statement checks the accelerationSettings of the specific body to see which acceleration is acting on the body. When a new acceleration model is made, a new case should be made for the new model. This case should have the same name as the name added to the AvailableAcceleration enum. Inside this case, the variable accelerationModelPointer should be assigned to the create acceleration model function, which was made before. When this is done, the new acceleration model should be incorporated into the acceleration framework of tudat.

## Adding a New State Derivative Model¶

Another type of model which allows the user to add extra features to it, is the state derivative model. These models are used by the dynamics simulator to determine how the state equations are solved. Currently there are several available, e.g. cowell state derivative model, encke state derivate model, and more, but if a special model is needed for a certain application, the user can add this to tudat by following this guide.

As before, this guide will start by looking at the framework of how the state derivative model is implemented. First, the state derivative model is chosen by the user from an enum called TranslationalPropagatorType, located in nBodyStateDerivativeModel.h. The model picked from this enum is then used as an input into the constructor of the TranslationalPropagatorSettings, where it is assigned to a specific member variable called: propagator_. The TranslationalPropagatorSettings is used as an input to the SingleArcDynamicsSimulator, where it will be used further. An example of this is shown below:

 // Create propagation settings.
std::shared_ptr< TranslationalStatePropagatorSettings< double > > propagatorSettings =
std::make_shared< TranslationalStatePropagatorSettings< double > >
( centralBodies, accelerationModelMap, bodiesToPropagate, systemInitialState,
terminationSettings, cowell, dependentVariablesToSave );

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


Inside SingleArcDynamicsSimulator, the propagator settings are passed to various functions, but the state derivative model is only needed in createStateDerivativeModels. This function returns a vector of SingleStateTypeDerivative. If a hybrid state derivative model is used, this function will fill the vector with state derivative models. If not, this vector will only contain one state derivative model, created by the function (equally) called createStateDerivativeModels. This function checks what kind of state is propagated (translational, rotational, etc.) and creates the specific state derivative model. For example, for the translational state it will call the createTranslationalStateDerivativeModel function. This function looks as follows:

 template< typename StateScalarType = double, typename TimeType = double >
std::shared_ptr< SingleStateTypeDerivative< StateScalarType, TimeType > >
createTranslationalStateDerivativeModel(
const std::shared_ptr< TranslationalStatePropagatorSettings< StateScalarType > >
translationPropagatorSettings,
const simulation_setup::NamedBodyMap& bodyMap,
const TimeType propagationStartTime )
{

// Create object for frame origin transformations.
std::shared_ptr< CentralBodyData< StateScalarType, TimeType > > centralBodyData =
createCentralBodyData< StateScalarType, TimeType >(
translationPropagatorSettings->centralBodies_,
translationPropagatorSettings->bodiesToIntegrate_,
bodyMap );

std::shared_ptr< SingleStateTypeDerivative< StateScalarType, TimeType > > stateDerivativeModel;

// Check propagator type and create corresponding state derivative object.
switch( translationPropagatorSettings->propagator_ )
{
case cowell:
{
stateDerivativeModel = std::make_shared<
NBodyCowellStateDerivative< StateScalarType, TimeType > >
( translationPropagatorSettings->getAccelerationsMap( ), centralBodyData,
translationPropagatorSettings->bodiesToIntegrate_ );
break;
}
...
...
...
}
default:
throw std::runtime_error(
"Error, did not recognize translational state propagation type: " +
std::to_string( translationPropagatorSettings->propagator_ ) );
}
return stateDerivativeModel;
}


This function is responsible for making the state derivate model based on the TranslationalPropagatorType variable decalred by the user, using a switch statement.

When a new state derivative model is implemented in tudat, two additions need to be made to the framework (besides making the model itself). The first part is to add the name of the model to the enum: TranslationalPropagatorType. This will allow the user to pick this model and use it as an input to the translational propagator settings. The second step is to add this model to the switch statement in createTranslationalStateDerivativeModel. A new case should be made for the name added to the enum before, and in it, the variable stateDerivativeModel should be assigned to the new model.

When building the new model, it is advised to use a state derivative model that is already available as an example for a place to start. The classes which contain these models are derived from the base class NBodyStateDerivative, and should contain a template for the TimeType that will be used. The NBodyStateDerivative class is again derived from another class called the SingleStateTypeDerivative. This class contains several pure virtual functions, which all should be added to the new model class in order for the new model to work. The specific names and input parameters of these functions can be found in the SingleStateTypeDerivative class. Once this is done, and the new model is implemented in the state derivative model framework, the new model should be available for the user.

Note

Don’t forget to put the include statement in createStateDerivativeModel.h if the new class is made in a seperate file.

## Adding a New Estimatable Parameter¶

The list of estimatable parameters already available in Tudat is presented in Setting Up Estimated Parameters. However, it is possible to add another parameter to this list of estimatable parameters if needed. This process requires to modify several files located in different directories and will be described in details based on the examples of the two parameters radiation_pressure_coefficient and rotation_pole_position.

First of all, the name of the new estimatable parameter has to be added to the list of the estimatable parameters available in Tudat, in the file estimatableParameter.h:

//! List of parameters that can be estimated by the orbit determination code
enum EstimatebleParametersEnum
{
arc_wise_initial_body_state,
initial_body_state,
initial_rotational_body_state,
gravitational_parameter,
constant_drag_coefficient,
spherical_harmonics_cosine_coefficient_block,
spherical_harmonics_sine_coefficient_block,
constant_rotation_rate,
rotation_pole_position,
...
...
}


In addition to the file estimatableParameters.h, the file estimatableParameters.cpp also has to be modified. In particular, a short description of each estimatable parameter has to be provided in the function getParameterTypeString, as it is done in the following for the parameters radiation_pressure_coefficient and rotation_pole_position:

std::string getParameterTypeString( const EstimatebleParametersEnum parameterType )
{
...
parameterDescription = "radiation pressure coefficient ";
break;
...
case rotation_pole_position:
parameterDescription = "pole position ";
break;
...
}


The type of the new estimatable parameter must then be specified within the function isDoubleParameter (still inside the file estimatableParameters.cpp). An estimatable parameter can either be a double or a Eigen::VectorXd. Regarding the two examples which are considered here, the parameter radiation_pressure_coefficient is a double while the parameter rotation_pole_position is a Eigen::VectorXd, whose first element is the right ascension of the rotation pole and the second one its declination.

bool isDoubleParameter( const EstimatebleParametersEnum parameterType )
{
...
isDoubleParameter = true;
break;
...
case rotation_pole_position:
isDoubleParameter = false;
break;
...
}


Finally, depending on which estimatable parameter is to be added to the list of available parameters, the functions isParameterDynamicalPropertyInitialState, isParameterRotationMatrixProperty, isParameterObservationLinkProperty and isParameterTidalProperty might also need to be modified to take this new parameter into account (again in estimatableParameters.cpp). They usually return a boolean which is set to false as default value, but a specific case has to added for the new parameter if it is related to either the initial state, a rotation matrix, an observation link or any tidal property. The parameter radiation_pressure_coefficient is related to none of them but the parameter rotation_pole_position is linked to a rotation matrix so that the function isParameterRotationMatrixProperty includes a case switching the boolean to true for this parameter.

bool isParameterRotationMatrixProperty( const EstimatebleParametersEnum parameterType )
{
bool flag;
switch( parameterType )
{
...
case rotation_pole_position:
flag = true;
break;
}
return flag;
}


Once the new estimatable parameter has been defined and characterised, a corresponding class is to be created to fully describe this parameter. This is done in an additional file which has to be created in the directory:

…/tudatBundle/tudat/Tudat/Astrodynamics/OrbitDetermination/EstimatableParameters/

Each estimatable parameter class is defined in a similar way, usually from the base class EstimatableParameter and includes the definition of three functions getParameterValues, setParameterValue and getParameterSize. Considering the parameter radiation_pressure_coefficient, the definition of its associated class RadiationPressureCoefficient is done as follows (in the file radiationPressureCoefficient.h).

class RadiationPressureCoefficient: public EstimatableParameter< double >
{

public:
//! Constructor.
{ }

//! Destructor.

//! Function to get the current value of the radiation pressure coefficient that is to be estimated.
double getParameterValue( )
{
}

//! Function to reset the value of the radiation pressure coefficient that is to be estimated.
void setParameterValue( double parameterValue )
{
}

//! Function to retrieve the size of the parameter.
int getParameterSize( ){ return 1; }

protected:

private:

//! Object containing the radiation pressure coefficient to be estimated.
};


For the parameter rotation_pole_position, the class ConstantRotationalOrientation is created in a very similar way, in a file constantRotationalOrientation.h, keeping in mind that this parameter is a Eigen::VectorXd whose size is 2 and not a double as it is the case for the radiation pressure coefficient.

Then an interface object has to be created to estimate the new parameter. This is done in the file createEstimatableParameters.cpp, either with the function createDoubleParameterToEstimate if the estimatable parameter is of type double or with createVectorParameterToEstimate if it is a Eigen::VectorXd object. These functions take as input an EstimatableParameterSettings object from which the settings of the parameter to be estimated are checked to ensure that it is properly defined and thus make the estimation possible. In some cases, the base class EstimatableParameterSettings is not detailed enough to give access to all the required properties of the estimatable parameter and a specific class has to be defined. As an example, the parameter spherical_harmonics_cosine_coefficient_block requires the definition of the estimatable parameter settings class called SphericalHarmonicEstimatableParameterSettings (defined in the file estimatableParameterSettings.h). This class manages the different combinations of degrees and orders of the spherical harmonic coefficients that have to be estimated.

Regarding the verification of the estimatable parameter settings of radiation_pressure_coefficient, it is checked that only one single radiation pressure interface is defined before creating the parameter and linking it to this radiation pressure interface and to the propagated body.

std::shared_ptr< EstimatableParameter< double > > createDoubleParameterToEstimate(
const std::shared_ptr< EstimatableParameterSettings >& doubleParameterName,
const NamedBodyMap& bodyMap, const basic_astrodynamics::AccelerationMap& accelerationModelMap )
{
...
{
if( currentBody->getRadiationPressureInterfaces( ).size( ) == 0 )
{
std::string errorMessage = "Error, no radiation pressure interfaces found in body " +
currentBodyName + " when making Cr parameter.";
throw std::runtime_error( errorMessage );
}
else if( currentBody->getRadiationPressureInterfaces( ).size( ) > 1 )
{
std::string errorMessage = "Error, multiple radiation pressure interfaces found in body " +
currentBodyName + " when making Cr parameter.";
throw std::runtime_error( errorMessage );
}
else
{
currentBodyName );
}
break;
}
...
}


Concerning the parameter rotation_pole_position, it must be verified that the rotation model is a simple rotational ephemeris for which the position of the rotation pole is indeed defined before creating the estimatable parameter.

std::shared_ptr< EstimatableParameter< Eigen::VectorXd > > createVectorParameterToEstimate(
const std::shared_ptr< EstimatableParameterSettings >& vectorParameterName,
const NamedBodyMap& bodyMap, const basic_astrodynamics::AccelerationMap& accelerationModelMap )
{
...
case rotation_pole_position:
if( std::dynamic_pointer_cast< SimpleRotationalEphemeris >( currentBody->getRotationalEphemeris( ) ) == nullptr )
{
std::string errorMessage = "Warning, no simple rotational ephemeris present in body " + currentBodyName +
" when making constant rotation orientation parameter";
throw std::runtime_error( errorMessage );
}
else
{
vectorParameterToEstimate = std::make_shared< ConstantRotationalOrientation >
( std::dynamic_pointer_cast< ephemerides::SimpleRotationalEphemeris >
( currentBody->getRotationalEphemeris( ) ), currentBodyName );

}
break;
...
}


To allow the parameter estimation to be conducted, partials with respect to this parameter have to be implemented. They can either be cartesian state or acceleration partials. In the former case, the partials of the cartesian state of the propagated body with respect to the estimatable parameter have to be implemented in the file createCartesianStatePartial.cpp. These partials are defined within the functions createCartesianStatePartialsWrtParameter (two functions exist with two different input types, depending on the type of the parameter (double or Eigen::VectorXd) that is to be considered). Three cases have to be distinguished here. First, if the estimatable parameter has been identified as being a property of a rotation matrix, then the function createCartesianStatePartialsWrtParameter calls another function named createRotationMatrixPartialsWrtParameter and defined in createCartesianStatePartial.cpp too (again two functions with the same name exist for the two types of estimatable parameters). A specific case has to be added within this function for each parameter which is related to a rotation matrix. For the parameter rotation_pole_position, the following lines of code have been added to first check that the rotation model is consistent with the estimatable parameter (here that is a simple rotational model) and then to call a function that returns the required rotation matrix partials for this particular model.

 std::shared_ptr< RotationMatrixPartial > createRotationMatrixPartialsWrtParameter(
const simulation_setup::NamedBodyMap& bodyMap,
const std::shared_ptr< estimatable_parameters::EstimatableParameter< Eigen::VectorXd > > parameterToEstimate )

{
...
case estimatable_parameters::rotation_pole_position:
if( std::dynamic_pointer_cast< ephemerides::SimpleRotationalEphemeris >(
currentBody->getRotationalEphemeris( ) ) == nullptr )
{
std::string errorMessage = "Warning, body's rotation model is not simple when making position w.r.t. pole position partial";
throw std::runtime_error( errorMessage );
}
// Create rotation matrix partial object
rotationMatrixPartial = std::make_shared< RotationMatrixPartialWrtPoleOrientation >(
std::dynamic_pointer_cast< SimpleRotationalEphemeris>( currentBody->getRotationalEphemeris( ) ) );
break;

...
}


The class RotationMatrixPartialWrtPoleOrientation used in the code above to define the rotationMatrixPartial variable is created in the file rotationMatrixPartial.h and must contain two internal functions calculatePartialOfRotationMatrixToBaseFrameWrtParameter and calculatePartialOfRotationMatrixDerivativeToBaseFrameWrtParameter. These functions return the rotation matrix and rotation matrix derivative partials respectively, with respect to the estimatable parameter. Specific functions to calculate these partials have to be added to the file rotationMatrixPartial.cpp. Regarding the parameter rotation_pole_position, these functions are calculatePartialOfRotationMatrixFromLocalFrameWrtPoleOrientation and calculatePartialOfRotationMatrixDerivativeFromLocalFrameWrtPoleOrientation, respectively.

So far, we have only considered the case where the estimatable parameter is related to a rotation matrix. However, when it is not the case, the function createCartesianStatePartialsWrtParameter is to be modified in a different way. A specific case has to be created for each parameter that is not a rotation matrix property. If the parameter has a direct impact on the cartesian state of the propagated body (eg ground_station_position), the partials of the cartesian state with respect to the parameters must be directly returned by the function.

A last case arises when the estimatable parameter neither is a rotation matrix property nor has a direct influence on the cartesian state of the propagated body but is involved in one of the acceleration models. The function createCartesianStatePartialsWrtParameter does not return any partial in that case. The estimation of such a parameter then requires the use of acceleration partials with respect to this parameter. The acceleration partials are implemented in the files of the directory:

…/tudatBundle/tudat/Tudat/Astrodynamics/OrbitDetermination/AccelerationPartials/

For any estimatable parameter related to an acceleration model, one must define a function getParameterPartialFunction which is a method of its associated acceleration partial clas. This function provides the partial of the acceleration model with respect to the estimatable parameter. Looking at the parameter radiation_pressure_coefficient in particular, the acceleration partials are defined from the class CannonBallRadiationPressurePartial which is itself derived from the base class AccelerationPartial). This acceleration partials class is defined in the files radiationPressureAccelerationPartial.h and radiationPressureAccelerationPartial.cpp. Only the radiation pressure coefficient can be estimated for this acceleration model so that the function getParameterPartialFunction only contains a single specific case dedicated to this parameter.

//! Function for setting up and retrieving a function returning a partial w.r.t. a double parameter.
std::pair< std::function< void( Eigen::MatrixXd& ) >, int > CannonBallRadiationPressurePartial::getParameterPartialFunction(
std::shared_ptr< estimatable_parameters::EstimatableParameter< double > > parameter )
{
std::function< void( Eigen::MatrixXd& ) > partialFunction;
int numberOfRows = 0;

// Check if parameter dependency exists.
if( parameter->getParameterName( ).second.first == acceleratedBody_ )
{
switch( parameter->getParameterName( ).first )
{
// Set function returning partial w.r.t. radiation pressure coefficient.

numberOfRows = 1;

break;
default:
break;
}
}
return std::make_pair( partialFunction, numberOfRows );
}


The function wrtRadiationPressureCoefficient called in the piece of code above is also a method of the class CannonBallRadiationPressurePartial and is defined in the file radiationPressureAccelerationPartial.h. It directly returns the partial derivative of the radiation pressure acceleration model with respect to the radiation pressure coefficient, as follows:

void wrtRadiationPressureCoefficient( Eigen::MatrixXd& partial )
{