Custom models

Attention

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 propagation models:

  • Custom thrust magnitude custom_thrust_magnitude()

  • Custom torque custom_torque()

  • Custom acceleration custom_acceleration()

  • Custom mass rate custom_mass_rate()

  • Custom state custom_state() Note: this custom function has current time, and the current custom state, as input

  • Custom termination setting custom_termination() Note: this custom function takes no inputs, if the current time is needed, it can be retrieved from the current_time attribute of the SystemOfBodies

  • 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_time attribute of the SystemOfBodies

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 getAerodynamicAngles and 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_guidance function 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