diff --git a/propagation/coupled_translational_rotational_dynamics.py b/propagation/coupled_translational_rotational_dynamics.py index e9d34ee..7f8f2eb 100644 --- a/propagation/coupled_translational_rotational_dynamics.py +++ b/propagation/coupled_translational_rotational_dynamics.py @@ -811,7 +811,7 @@ def get_longitudinal_normal_mode_from_inertia_tensor(inertia_tensor: np.ndarray, """ As we just saw, Phobos does not rotate the way we would expect it to, due to the presence of excited normal modes. The fact that the normal mode is excited is a result of the initial rotational state, which is not consistent with rotational dynamics in which the normal modes are damped. We solve this problem by determining an initial state where these normal modes **are** damped, following the method used by (Rambaux et al. 2010). Specifically, we introduce a _virtual torque_ ("virtual" is a scientific word for "made up") that we know has a strong damping effect, and we propagate the dynamics forward in time. At one point, we will reach a state in which the normal modes have been eliminated from the dynamics (or 'damped'). Now, we remove this damping torque and we propagate the dynamics backward in time from this point onwards, until we reach our original initial time. This new rotational state that we have found, and will call the _damped initial state_ would ideally not excite the normal modes of Phobos when propagated forward again using our 'normal' (e.g. without the virtual torque) dynamical model. -In practice, the normal mode will be excited much less than initially, but some excitation still remains. This is further resolved by repeating the above procedure many times, with increasing 'damping times'. In Tudat, this approach is automated by the `get_zero_proper_mode_rotational_state` function [TODO: add API docs], where additional details of the algorithm are described. In short, a _damping time_ :math:`\tau_{d}` is selected, which is a measure of how quick the virtual torque will eliminate the normal modes, and it is used to compute the torque. The propagation is performed forwards in time with the virtual torque enabled, for :math:`10 \tau_{d}`, and backwards to the original time without it. The new initial state is used in another iteration with a new and larger damping time. The implementation of the virtual torque in this algorithm is very specific to bodies in spin-orbit resonance, which will (averaged over long time periods) rotate uniformly around their body-fixed :math:`z` axis, with a period equal to the orbital period. The damping torque will act to prevent the excitation of the normal modes while preserving that uniform rotation. The function thus requires as inputs: (1) the `bodies` object, (2) the `propagator_settings` object, (3) the mean rotational rate of the body around its Z axis, and (4) a list of dissipation times or damping times in order to compute the virtual torques for the different iteration. The above is implemented as follows: +In practice, the normal mode will be excited much less than initially, but some excitation still remains. This is further resolved by repeating the above procedure many times, with increasing 'damping times'. In Tudat, this approach is automated by the `get_damped_proper_mode_initial_rotational_state` function [TODO: add API docs], where additional details of the algorithm are described. In short, a _damping time_ :math:`\tau_{d}` is selected, which is a measure of how quick the virtual torque will eliminate the normal modes, and it is used to compute the torque. The propagation is performed forwards in time with the virtual torque enabled, for :math:`10 \tau_{d}`, and backwards to the original time without it. The new initial state is used in another iteration with a new and larger damping time. The implementation of the virtual torque in this algorithm is very specific to bodies in spin-orbit resonance, which will (averaged over long time periods) rotate uniformly around their body-fixed :math:`z` axis, with a period equal to the orbital period. The damping torque will act to prevent the excitation of the normal modes while preserving that uniform rotation. The function thus requires as inputs: (1) the `bodies` object, (2) the `propagator_settings` object, (3) the mean rotational rate of the body around its Z axis, and (4) a list of dissipation times or damping times in order to compute the virtual torques for the different iteration. The above is implemented as follows: """ phobos_mean_rotational_rate = 0.000228035245 # In rad/s diff --git a/propagation/thrust_satellite_engine.ipynb b/propagation/thrust_satellite_engine.ipynb index 0a57b16..8cecae0 100644 --- a/propagation/thrust_satellite_engine.ipynb +++ b/propagation/thrust_satellite_engine.ipynb @@ -348,7 +348,7 @@ "outputs": [], "execution_count": null, "source": [ - "simulation_start_date = DateTime(2035,7,28,14,24)\n", + "simulation_start_date = time_conversion.DateTime(2035,7,28,14,24)\n", "simulation_start_epoch = simulation_start_date.epoch()\n", "simulation_end_epoch = simulation_start_epoch + 344.0 * constants.JULIAN_DAY / 24.0\n", "\n", diff --git a/propagation/thrust_satellite_engine.py b/propagation/thrust_satellite_engine.py index 9cb1499..2042737 100644 --- a/propagation/thrust_satellite_engine.py +++ b/propagation/thrust_satellite_engine.py @@ -1,36 +1,37 @@ -#!/usr/bin/env python -# coding: utf-8 -# # Basic thrust model for satellite engines -# -# Copyright (c) 2010-2022, Delft University of Technology. All rights reserved. This file is part of the Tudat. Redistribution and use in source and binary forms, with or without modification, are permitted exclusively under the terms of the Modified BSD license. You should have received a copy of the license with this file. If not, please visit: http://tudat.tudelft.nl/LICENSE. +# Basic thrust model for satellite engines +""" -# ## Context -# -# This example demonstrates the application of a basic custom class to model the thrust of a satellite. -# -# Consider the orbit of JUICE with respect to Ganymede. The goal will be to increase the inclination of the orbit, while keeping the nodes constant. Theoretically speaking, the easiest way to do this is by one or more impulsive maneuvers at the nodes of the orbit. Naturally, in reality this thrust will not be impulsive, but the thrust will be applied over a finite duration. -# -# In this example, thrust is implemented to be provided when the true longitude of the spacecraft is within 2 degrees of one of the nodes. Note that parametrizing the thrust as a function of the true longitude is likely not an optimal choice, but it will suffice for this example. The maximum thrust magnitude $T_{max}$ is used when the thrust is on. Introducing such a discontinuity in the dynamical model is again likely not an optimal choice, but will suffice. -# In addition to the acceleration from the thrust, a basic model is set up, consisting of the following accelerations: \ -# · Spherical harmonic gravity accelerations from Ganymede and Jupiter. Their gravity fields will be expanded to D/O 2/2 and 4/0 respectively. \ -# · Third body forces will include those of the Sun, Saturn, Europa, Io and Callisto. -# -# The mass of the vehicle is propagated using a mass rate model consistent with the engine thrust used. -# -# How to set up dependent variables to save the mass and altitude of the vehicle over time is demonstrated. -# -# Finally, some post-processing is performed to analyze the retrieved results. +Copyright (c) 2010-2022, Delft University of Technology. All rights reserved. This file is part of the Tudat. Redistribution and use in source and binary forms, with or without modification, are permitted exclusively under the terms of the Modified BSD license. You should have received a copy of the license with this file. If not, please visit: http://tudat.tudelft.nl/LICENSE. +""" -# ## Import statements -# -# The required import statements are made here at the beginning of the code. -# -# Some standard modules are first loaded: `numpy` and `matplotlib.pyplot`. -# -# The different modules of `tudatpy` that will be used are then imported. We include here the loading of the standard spice kernels as well as the JUICE spice kernel. This is not a standard kernel included in SPICE, and hence an associated file will need to be provided as input, containing the relevant data. +## Context +""" + +This example demonstrates the application of a basic custom class to model the thrust of a satellite. + +Consider the orbit of JUICE with respect to Ganymede. The goal will be to increase the inclination of the orbit, while keeping the nodes constant. Theoretically speaking, the easiest way to do this is by one or more impulsive maneuvers at the nodes of the orbit. Naturally, in reality this thrust will not be impulsive, but the thrust will be applied over a finite duration. + +In this example, thrust is implemented to be provided when the true longitude of the spacecraft is within 2 degrees of one of the nodes. Note that parametrizing the thrust as a function of the true longitude is likely not an optimal choice, but it will suffice for this example. The maximum thrust magnitude $T_{max}$ is used when the thrust is on. Introducing such a discontinuity in the dynamical model is again likely not an optimal choice, but will suffice. +In addition to the acceleration from the thrust, a basic model is set up, consisting of the following accelerations: \ +· Spherical harmonic gravity accelerations from Ganymede and Jupiter. Their gravity fields will be expanded to D/O 2/2 and 4/0 respectively. \ +· Third body forces will include those of the Sun, Saturn, Europa, Io and Callisto. + +The mass of the vehicle is propagated using a mass rate model consistent with the engine thrust used. + +How to set up dependent variables to save the mass and altitude of the vehicle over time is demonstrated. + +Finally, some post-processing is performed to analyze the retrieved results. +""" -# In[1]: +## Import statements +""" +The required import statements are made here at the beginning of the code. + +Some standard modules are first loaded: `numpy` and `matplotlib.pyplot`. + +The different modules of `tudatpy` that will be used are then imported. We include here the loading of the standard spice kernels as well as the JUICE spice kernel. This is not a standard kernel included in SPICE, and hence an associated file will need to be provided as input, containing the relevant data. +""" import numpy as np from matplotlib import pyplot as plt @@ -49,16 +50,17 @@ spice.load_kernel("input/juice_mat_crema_5_1_150lb_v01.bsp") -# ## Environment setup -# -# Let’s create the environment for our simulation. This setup covers the creation of (celestial) bodies, vehicle(s), and environment interfaces. +## Environment setup +""" -# ### Create the bodies -# -# The required bodies will first be created. Note that this involves creation of a simple atmosphere for Ganymede as well as including the radiation pressure settings for JUICE. Please refer to the [Environment Models](https://docs.tudat.space/en/latest/_src_user_guide/state_propagation/environment_setup/environment_models.html) in the user guide for more details. +Let’s create the environment for our simulation. This setup covers the creation of (celestial) bodies, vehicle(s), and environment interfaces. +""" -# In[3]: +### Create the bodies +""" +The required bodies will first be created. Note that this involves creation of a simple atmosphere for Ganymede as well as including the radiation pressure settings for JUICE. Please refer to the [Environment Models](https://docs.tudat.space/en/latest/_src_user_guide/state_propagation/environment_setup/environment_models.html) in the user guide for more details. +""" # Create settings for celestial bodies bodies_to_create = ["Ganymede", "Sun", "Jupiter", "Saturn", "Europa", "Io", "Callisto"] @@ -86,16 +88,15 @@ bodies = environment_setup.create_system_of_bodies(body_settings) -# ### Define the thrust guidance settings -# -# Let's now define the aspects of the body related to thrust: its orientation and its engine. This is done by defining a `ThrustGuidance` class that will contain functions that will calculate the thrust magnitude and thrust direction at any point in time. -# -# The class takes as input `maximum_thrust`, which defines what the constant value of the thrust will be when turned on in Newtons; `true_longitude_threshold`, which defines the range in degrees within which the thrust will be turned on, as measured from either of the nodes; `bodies` is the list of bodies created of type `SystemOfBodies`. Note that, in this class the thrust direction and magnitude are computed entirely independently (despite the fact that some computations are then done twice) to reduce the complexity of the example. -# -# Then, the `rotation_model_settings` and `thrust_magnitude_settings` are set up, which complete the thrust model. +### Define the thrust guidance settings +""" + +Let's now define the aspects of the body related to thrust: its orientation and its engine. This is done by defining a `ThrustGuidance` class that will contain functions that will calculate the thrust magnitude and thrust direction at any point in time. -# In[4]: +The class takes as input `maximum_thrust`, which defines what the constant value of the thrust will be when turned on in Newtons; `true_longitude_threshold`, which defines the range in degrees within which the thrust will be turned on, as measured from either of the nodes; `bodies` is the list of bodies created of type `SystemOfBodies`. Note that, in this class the thrust direction and magnitude are computed entirely independently (despite the fact that some computations are then done twice) to reduce the complexity of the example. +Then, the `rotation_model_settings` and `thrust_magnitude_settings` are set up, which complete the thrust model. +""" class ThrustGuidance: @@ -174,9 +175,6 @@ def compute_thrust_magnitude(self, current_time: float): return 0.0 -# In[5]: - - # Create thrust guidance object (e.g. object that calculates direction/magnitude of thrust) thrust_magnitude = 1 # N true_longitude_threshold = 2 # degrees @@ -200,15 +198,14 @@ def compute_thrust_magnitude(self, current_time: float): environment_setup.add_rotation_model(bodies, "JUICE", rotation_model_settings) -# ## Propagation setup -# -# Now that the environment is created, the propagation setup is defined. -# -# First, the bodies to be propagated and the central bodies will be defined. -# Central bodies are the bodies with respect to which the state of the respective propagated bodies is defined. +## Propagation setup +""" -# In[6]: +Now that the environment is created, the propagation setup is defined. +First, the bodies to be propagated and the central bodies will be defined. +Central bodies are the bodies with respect to which the state of the respective propagated bodies is defined. +""" # Define bodies that are propagated bodies_to_propagate = ["JUICE"] @@ -217,16 +214,15 @@ def compute_thrust_magnitude(self, current_time: float): central_bodies = ["Ganymede"] -# ### Create the acceleration model -# -# First off, the acceleration settings that act on the JUICE are to be defined, which consists of the accelerations previously mentioned. -# -# The defined acceleration settings are then applied to JUICE in a dictionary. -# -# Finally, this dictionary is input to the propagation setup to create the acceleration models. +### Create the acceleration model +""" + +First off, the acceleration settings that act on the JUICE are to be defined, which consists of the accelerations previously mentioned. -# In[7]: +The defined acceleration settings are then applied to JUICE in a dictionary. +Finally, this dictionary is input to the propagation setup to create the acceleration models. +""" # Define accelerations acting on vehicle. @@ -273,13 +269,13 @@ def compute_thrust_magnitude(self, current_time: float): bodies, acceleration_settings, bodies_to_propagate, central_bodies) -# ### Define the initial state -# -# The initial state of JUICE that will be propagated is now defined. -# -# In this example, the initial state is retrieved from the JUICE spice kernel for a certain epoch, as defined at the start of the script. The final epoch is also given, making the propagation last a little over two weeks. +### Define the initial state +""" -# In[23]: +The initial state of JUICE that will be propagated is now defined. + +In this example, the initial state is retrieved from the JUICE spice kernel for a certain epoch, as defined at the start of the script. The final epoch is also given, making the propagation last a little over two weeks. +""" simulation_start_date = time_conversion.DateTime(2035,7,28,14,24) simulation_start_epoch = simulation_start_date.epoch() @@ -294,16 +290,15 @@ def compute_thrust_magnitude(self, current_time: float): ephemeris_time=simulation_start_epoch) -# ### Define dependent variables to save -# -# In this example, we are interested in saving not only the propagated state of the satellite over time, but also a set of so-called dependent variables that are to be computed (or extracted and saved) at each integration step. -# -# [This page](https://tudatpy.readthedocs.io/en/latest/dependent_variable.html) of the tudatpy API website provides a detailled explanation of all the dependent variables that are available. -# -# We will save the keplerian state of JUICE with respect to Ganymede, to verify whether we really do raise the inclincation of the orbit. We will also save the mass of JUICE over time to ensure our mass propagation is also correctly performed. Lastly, we save norm of the thrust acceleration to verify whether the thrust we have set is indeed equal to 1N. +### Define dependent variables to save +""" -# In[24]: +In this example, we are interested in saving not only the propagated state of the satellite over time, but also a set of so-called dependent variables that are to be computed (or extracted and saved) at each integration step. +[This page](https://tudatpy.readthedocs.io/en/latest/dependent_variable.html) of the tudatpy API website provides a detailled explanation of all the dependent variables that are available. + +We will save the keplerian state of JUICE with respect to Ganymede, to verify whether we really do raise the inclincation of the orbit. We will also save the mass of JUICE over time to ensure our mass propagation is also correctly performed. Lastly, we save norm of the thrust acceleration to verify whether the thrust we have set is indeed equal to 1N. +""" # Define required outputs dependent_variables_to_save = [ @@ -315,27 +310,25 @@ def compute_thrust_magnitude(self, current_time: float): ] -# ### Create the termination settings -# -# Let's now define a set of termination settings. In this setup, we only define a single termination setting: -# - Stop when JUICE reaches the specified end epoch (after a little over 14 days). -# -# In a more realistic scenario, you may want to set multiple termination settings; next to terminating on time you can also choose to terminate at a certain mass (e.g. when you run out of propellant) or a certain altitude. For simplicity, we simply let the propagation run its full course and we will analyse the results after. +### Create the termination settings +""" -# In[25]: +Let's now define a set of termination settings. In this setup, we only define a single termination setting: +- Stop when JUICE reaches the specified end epoch (after a little over 14 days). +In a more realistic scenario, you may want to set multiple termination settings; next to terminating on time you can also choose to terminate at a certain mass (e.g. when you run out of propellant) or a certain altitude. For simplicity, we simply let the propagation run its full course and we will analyse the results after. +""" termination_settings = propagation_setup.propagator.time_termination(simulation_end_epoch) -# ### Create the integrator settings -# -# The last step before starting the simulation is to set up the integrator that will be used. -# -# In this case, an RK4 fixed step size integrator is used, with a given step size of 10 s. This is one of the most simple and reliable integrators available, and usually works good enough for a wide number of applications, provided there are relatively stable dynamics. +### Create the integrator settings +""" -# In[26]: +The last step before starting the simulation is to set up the integrator that will be used. +In this case, an RK4 fixed step size integrator is used, with a given step size of 10 s. This is one of the most simple and reliable integrators available, and usually works good enough for a wide number of applications, provided there are relatively stable dynamics. +""" # Create numerical integrator settings. fixed_step_size = 10.0 @@ -344,25 +337,27 @@ def compute_thrust_magnitude(self, current_time: float): ) -# ### Create the propagator settings -# -# The propagator settings are now defined. -# -# As usual, translational propagation settings are defined, using a Cowell propagator. -# -# Next, mass propagation settings are also defined. The mass rate from the vehicle is set up to be consistent with the thrust used. The vehicle then loses mass as it loses propellant. -# -# Finally, a multitype propagator is defined, encompassing both the translational and the mass propagators. -# -# ## Coupled dynamics -# If you have used Tudat before, you are most probably familiar with what _translational propagators_ are. Possibly, you are also familiar with combined translational-mass propagations. These are just an example of a **multi-type propagation**. The way Tudat deals with these multi-type propagations is by creating the appropriate "single-type" propagation settings for each type of dynamics, and then putting them all together at then end in the _multi-type propagator settings_. Thus, we will follow the same process here. For more details, see [multi-type propagation documentation](https://docs.tudat.space/en/latest/_src_user_guide/state_propagation/propagation_setup/multi_type.html). -# -# As you will see in the code below - and can be deduced comparing the APIs for the [translational](https://py.api.tudat.space/en/latest/propagator.html#tudatpy.numerical_simulation.propagation_setup.propagator.translational), [mass](https://py.api.tudat.space/en/latest/propagator.html#tudatpy.numerical_simulation.propagation_setup.propagator.mass) and [multi-type](https://py.api.tudat.space/en/latest/propagator.html#tudatpy.numerical_simulation.propagation_setup.propagator.multitype) propagators - some of the inputs (namely the integrator settings, the initial time, the termination settings and the output variables) are identical between all three propagators - the two single-type and the one multi-type. In these overlaps, tudat will only read the "top level" arguments, i.e. those passed to the multi-type propagator and will ignore the rest. This means that these inputs can be left empty (`0`, `NaN` or `None`) for the single-type propagators. However, it is good practice to be self-consistent and pass the same inputs to all propagators. This facilitates the use of the single-type propagators for the simulation of only one type of dynamics while being consistent with the inputs of the multi-type simulation. -# -# Below, we will begin by creating these common inputs and will then move on to the propagator-specific inputs. +### Create the propagator settings +""" + +The propagator settings are now defined. + +As usual, translational propagation settings are defined, using a Cowell propagator. + +Next, mass propagation settings are also defined. The mass rate from the vehicle is set up to be consistent with the thrust used. The vehicle then loses mass as it loses propellant. + +Finally, a multitype propagator is defined, encompassing both the translational and the mass propagators. + +""" + +## Coupled dynamics +""" +If you have used Tudat before, you are most probably familiar with what _translational propagators_ are. Possibly, you are also familiar with combined translational-mass propagations. These are just an example of a **multi-type propagation**. The way Tudat deals with these multi-type propagations is by creating the appropriate "single-type" propagation settings for each type of dynamics, and then putting them all together at then end in the _multi-type propagator settings_. Thus, we will follow the same process here. For more details, see [multi-type propagation documentation](https://docs.tudat.space/en/latest/_src_user_guide/state_propagation/propagation_setup/multi_type.html). -# In[27]: +As you will see in the code below - and can be deduced comparing the APIs for the [translational](https://py.api.tudat.space/en/latest/propagator.html#tudatpy.numerical_simulation.propagation_setup.propagator.translational), [mass](https://py.api.tudat.space/en/latest/propagator.html#tudatpy.numerical_simulation.propagation_setup.propagator.mass) and [multi-type](https://py.api.tudat.space/en/latest/propagator.html#tudatpy.numerical_simulation.propagation_setup.propagator.multitype) propagators - some of the inputs (namely the integrator settings, the initial time, the termination settings and the output variables) are identical between all three propagators - the two single-type and the one multi-type. In these overlaps, tudat will only read the "top level" arguments, i.e. those passed to the multi-type propagator and will ignore the rest. This means that these inputs can be left empty (`0`, `NaN` or `None`) for the single-type propagators. However, it is good practice to be self-consistent and pass the same inputs to all propagators. This facilitates the use of the single-type propagators for the simulation of only one type of dynamics while being consistent with the inputs of the multi-type simulation. +Below, we will begin by creating these common inputs and will then move on to the propagator-specific inputs. +""" # Create propagation settings. translational_propagator_settings = propagation_setup.propagator.translational( @@ -405,24 +400,23 @@ def compute_thrust_magnitude(self, current_time: float): propagator_settings.print_settings.print_initial_and_final_conditions = True -# ## Propagate the orbit -# -# The orbit of JUICE around Ganymede is now ready to be propagated. -# -# This is done by calling the `create_dynamics_simulator()` function of the `numerical_simulation module`. -# This function requires the `system_of_bodies` and `propagator_settings` that have all been defined earlier. -# -# After this, the history of the propagated state over time, containing both the position and velocity history, is extracted. -# This history, taking the form of a dictionary, is then converted to an array containing 7 columns: -# - Column 0: Time history, in seconds since J2000. -# - Columns 1 to 3: Position history, in meters, in the frame that was specified in the `body_settings`. -# - Columns 4 to 6: Velocity history, in meters per second, in the frame that was specified in the `body_settings`. -# -# The same is done with the dependent variable history. The column indexes corresponding to a given dependent variable in the `dependent_variables_to_save` variable are printed when the simulation is run, when `create_dynamics_simulator()` is called. -# Do pay attention that converting to an `ndarray` using the `result2array()` utility will shift these indexes because the first column (index 0) will then be the times. +## Propagate the orbit +""" -# In[28]: +The orbit of JUICE around Ganymede is now ready to be propagated. +This is done by calling the `create_dynamics_simulator()` function of the `numerical_simulation module`. +This function requires the `system_of_bodies` and `propagator_settings` that have all been defined earlier. + +After this, the history of the propagated state over time, containing both the position and velocity history, is extracted. +This history, taking the form of a dictionary, is then converted to an array containing 7 columns: +- Column 0: Time history, in seconds since J2000. +- Columns 1 to 3: Position history, in meters, in the frame that was specified in the `body_settings`. +- Columns 4 to 6: Velocity history, in meters per second, in the frame that was specified in the `body_settings`. + +The same is done with the dependent variable history. The column indexes corresponding to a given dependent variable in the `dependent_variables_to_save` variable are printed when the simulation is run, when `create_dynamics_simulator()` is called. +Do pay attention that converting to an `ndarray` using the `result2array()` utility will shift these indexes because the first column (index 0) will then be the times. +""" # Create simulation object and propagate dynamics. dynamics_simulator = numerical_simulation.create_dynamics_simulator( @@ -436,14 +430,13 @@ def compute_thrust_magnitude(self, current_time: float): dependent_variable_history = propagation_results.dependent_variable_history -# ## Let's look at plots -# -# We are now ready to do some post-processing, and look at our results! As mentioned before, we are interested in validating whether our maneuver was performed successfully, and whether the mass propagation has accurately kept track of the mass of JUICE over time. -# -# As such, we'll plot the inclination, the mass of JUICE over time and add several plots to gain some confidence in the correct implementation of the thrust direction and magnitude, such as a plot of the magnitude as a function of the true longitude. +## Let's look at plots +""" -# In[29]: +We are now ready to do some post-processing, and look at our results! As mentioned before, we are interested in validating whether our maneuver was performed successfully, and whether the mass propagation has accurately kept track of the mass of JUICE over time. +As such, we'll plot the inclination, the mass of JUICE over time and add several plots to gain some confidence in the correct implementation of the thrust direction and magnitude, such as a plot of the magnitude as a function of the true longitude. +""" ## PLOTS