This page is only compatible for Tudatpy version >= 0.7. See thrust and rotation refactor for more information on this page for older versions.
When working with a very specific model or application, it often happens that the model you want to use is not implemented in Tudat. Two important examples are:
If this is the case, we have the option for users to define ‘custom’ models for various models in both the environment and propagation setup modules. The use of these custom settings requires the user to define their own function for the specific model, as is shown on this page with a number of examples. Below, you can find a list of the currently supported custom environment models in Tudat. In most cases, the input to the custom function is the current time only. In case the user wants to define other state/environment dependencies, these can be implemented using a user-defined custom class, as shown below. Custom functions that have additional dependencies are noted explicitly below. The output of the custom function will differ per model type (e.g. density custom model: float output; state custom model: vector output), and are specified in the API entries that are linked to.
Custom environment models:
Custom aerodynamic coefficients (with suppported dependencies from
custom_aerodynamic_force_and_moment_coefficients()(only force, and force and moment coefficients, respectively)
Custom atmospheric density model (function of time only)
Custom atmospheric density model (function of position and time)
custom_four_dimensional_constant_temperature()Note: this custom function has current time, altitude, latitude and longitude as input.
Custom wind model:
Custom body orientation
Custom inertial direction (typical for basic thrust guidance)
Custom propagation models:
Custom thrust magnitude
Custom mass rate
custom_state()Note: this custom function has current time, and the current custom state, as input
Custom dependent variable
custom_dependent_variable()Note: this custom function takes no inputs, if the current time is needed, it can be retrieved from the
current_timeattribute of the
In each case, the user is required to define their own function, with a predefined set of inputs and outputs, which are different for each specific environment model (see API documentation links above). We can distinguish (roughly) three different ways in which to provide such custom functions to Tudat:
As a Python free function, which is responsible for calculating a single custom model
As a member function of a Python class, which is responsible for calculating a single custom model
As a member function of a Python class, which is responsible for calculating a multiple custom model
The latter choice permits complex guidance algorithms to be implemented, in which (for instance) the algorithm for the control surface deflection, thrust magnitude and body-fixed thrust direction (thrust vector control) are calculated in a coupled manner.
For most environment models, the custom model can be fully defined by a standalone function , and can be fully defined by a free function (which may of course call other functions from tudat, other packages or your own code. A custom function will, in numerous cases, have to depend on properties of bodies (states, altitude, density, etc.). How to access such properties of the environment during the propagation is described on the page on Interacting with the environment during propagation.
Custom model, free function
Here, we show an example of an ephemeris model that will be both faster, and less accurate, than the models implemented in Tudat. This may be advantageous for preliminary simulation. The model is very simple: a perfectly circular orbit in the \(xy-\) frame.
def neptune_state_function( current_time ): if( current_time == current_time ): # Define constants for Neptune ephemeris orbital_radius = 4.5E12 sun_gravitational_parameter = 1.32712440042E20 anomaly_at_j2000 = np.deg2rad( 304.88003 ) # Compute orbitall velocity and orbital period orbital_velocity = math.sqrt( sun_gravitational_parameter / orbital_radius ) orbital_period = 2.0 * math.pi * math.sqrt( orbital_radius ** 3 / sun_gravitational_parameter ) # Compute current angular position along orbit anomaly_at_epoch = anomaly_at_j2000 + 2.0 * math.pi * ( current_time / orbital_period ) # Compute and return Neptune state return np.array([orbital_radius * math.cos(anomaly_at_epoch)],[orbital_radius * math.cos(anomaly_at_epoch)],[0.0] ) ... # Retrieve custom state function, and set as ephemeris function w.r.t. Sun, with axes along J2000 custom_state_function = neptune_state_function body_settings.get( "Neptune" ).ephemeris_settings = environment_setup.ephemeris.custom( custom_state_function, 'Sun', 'ECLIPJ2000' )
In the above example, the user-define function
neptune_state_function is provided when creating the custom ephemeris settings. The only requirement on this custom function is that it takes a single float as argument (representing time since J2000), and returns a 6-dimensional vector (representing the Cartesian state in the frame specified), as can be seen in the
custom() API entry.
Custom model (single), class member
Below, a skeleton is given for a custom class for the calculation of the thrust magnitude.
class SimpleCustomGuidanceModel: def __init__(self, bodies: environment.SystemOfBodies): # Extract the STS and Earth bodies self.vehicle = bodies.get_body("Vehicle") self.earth = bodies.get_body("Earth") def getThrustMagnitude(self, current_time: float): if( current_time == current_time ): # Update the class to the current time self.thrust_magnitude = ... # Return angles calculated by update function return self.thrust_magnitude
This setup allows the guidance model to direcly access any of the properties of the bodies named ‘Earth’ and ‘Vehicle’, which were set as class attributes of the
SimpleCustomGuidanceModel class (note: the inputs to the constructor, and the manner in which they are used is entirely up to the user, here we give just one example).
The custom thrust magnitude model can then be used as follows to define the thrust magnitude that is to be exerted by an engine:
# Create bodies (same as in any other simulation) bodies = ... # Create guidance object guidance_model = SimpleCustomGuidanceModel( bodies ) # Extract guidance function thrust_magnitude_function = guidance_model.get_thrust_magnitude # Create thrust settings from custom model thrust_magnitude_settings = environment_setup.thrust.custom_thrust_magnitude( thrust_magnitude_function, specific_impulse = 300.0 ) # Create engine model, and add to Vehicle environment_setup.add_engine_model( "Vehicle", "MainEngine", thrust_magnitude_settings, system_of_bodies )
Custom model (multiple), class member
Below, a skeleton is given for a custom class for the calculation of both thrust magnitude and aerodynamic angles.
class CustomGuidanceModel: def __init__(self, bodies: environment.SystemOfBodies): # Extract the STS and Earth bodies self.vehicle = bodies.get_body("Vehicle") self.earth = bodies.get_body("Earth") self.current_time = float("NaN") def get_aerodynamic_angles(self, current_time: float): # Update the class to the current time self.update_guidance( current_time ) # Return angles calculated by update function return np.array([self.angle_of_attack, 0.0, self.bank_angle]) def get_thrust_magnitude(self, current_time: float): # Update the class to the current time self.update_guidance( current_time ) # Return angles calculated by update function return self.thrust_magnitude def update_guidance(self, current_time: float): if( math.isnan( current_time ) ): # Set the model's current time to NaN, indicating that it needs to be updated self.current_time = float("NaN") elif( current_time != self.current_time ): # Calculate current body orientation through angle of attack and bank angle self.angle_of_attack = ... self.bank_angle = ... # Calculate current thrust magnitude self.thrust_magnitude = ... # Set the model's current time, indicating that it has been updated self.current_time = current_time
Here, we see a different setup compared to the previous case. There is a single function, named
update_guidance that calculates both the thrust magnitude and the aerodynamic angles. This allows the calculation of the two models to be coupled, which is required for many more advanced applications. The two functions
getThrustMagnitude are then linked to the environment as follows:
# Create bodies (same as in any other simulation) bodies = ... # Create guidance object guidance_model = SimpleCustomGuidanceModel( bodies ) # Extract guidance function thrust_magnitude_function = guidance_model.get_thrust_magnitude aerodynamic_angle_function = guidance_model.get_aerodynamic_angles # Create thrust settings from custom model, create engine model and add it to vehicle thrust_magnitude_settings = environment_setup.thrust.custom_thrust_magnitude( thrust_magnitude_function, specific_impulse = 300.0 ) environment_setup.add_engine_model( "Vehicle", "MainEngine", thrust_magnitude_settings, system_of_bodies ) # Create angle-based rotation model settings for vehicle, create rotation model, and add it to vehicle rotation_model_settings = environment_setup.rotation_model.aerodynamic_angle_based( "Earth", "J2000", "Vehicle_fixed", aerodynamic_angle_function) environment_setup.add_rotation_model( bodies, "Vehicle", rotation_model_settings )
In setting up the custom guidance class, we now need to take care of one crucial point: even though data is retrieved from the object twice per function evaluation of the state derivative, the calculation should only be done once. Since it is often difficult to predict which of the custom functions will be called first, we use a different setup: defining a
current_time member variable, and letting the code check whether an update needs to be done. This is achieved as follows:
After the guidance function is evaluated, the class member time is set to the input time, and the guidance is not evaluated a second time during the same state derivative function evaluation
At the very start of a state derivative function evaluation, the
update_guidancefunction is called with a NaN input (done by each custom function) signalling that a new function evaluation has started, and the class needs to recompute the guidance. This is done to support integrators such as the RK4 integrator, where two successive state derivatives are evaluated using the same time, but different states
If the current time of the class is NaN, the guidance is by definition recomputed when called