diff --git a/CHANGELOG.md b/CHANGELOG.md
index 87e3c50a4..4d2509adf 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -32,11 +32,15 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
 
 ### Added
 
+- ENH: Exponential backoff decorator (fix #449) [#588](https://github.com/RocketPy-Team/RocketPy/pull/588)
 - ENH: Add new stability margin properties to Flight class [#572](https://github.com/RocketPy-Team/RocketPy/pull/572)
 - ENH: adds `Function.remove_outliers` method [#554](https://github.com/RocketPy-Team/RocketPy/pull/554)
 
 ### Changed
 
+- DEP: delete deprecated rocketpy.tools.cached_property [#587](https://github.com/RocketPy-Team/RocketPy/pull/587) 
+- MNT: Modularize Rocket Draw [#580](https://github.com/RocketPy-Team/RocketPy/pull/580)
+- DOC: Improvements of Environment docstring phrasing [#565](https://github.com/RocketPy-Team/RocketPy/pull/565)
 - MNT: Refactor flight prints module [#579](https://github.com/RocketPy-Team/RocketPy/pull/579)
 - DOC: Convert CompareFlights example notebooks to .rst files [#576](https://github.com/RocketPy-Team/RocketPy/pull/576)
 - MNT: Refactor inertia calculations using parallel axis theorem [#573] (https://github.com/RocketPy-Team/RocketPy/pull/573)
@@ -44,6 +48,7 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
 
 ### Fixed
 
+- BUG: Swap rocket.total_mass.differentiate for motor.total_mass_flow rate [#585](https://github.com/RocketPy-Team/RocketPy/pull/585)
 - BUG: export_eng 'Motor' method would not work for liquid motors. [#559](https://github.com/RocketPy-Team/RocketPy/pull/559) 
 
 ## [v1.2.2] - 2024-03-22
diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py
index 546d446e8..5a845a173 100644
--- a/rocketpy/environment/environment.py
+++ b/rocketpy/environment/environment.py
@@ -3,7 +3,7 @@
 import re
 import warnings
 from collections import namedtuple
-from datetime import datetime, timedelta
+from datetime import datetime, timedelta, timezone
 
 import numpy as np
 import numpy.ma as ma
@@ -13,6 +13,7 @@
 from ..mathutils.function import Function, funcify_method
 from ..plots.environment_plots import _EnvironmentPlots
 from ..prints.environment_prints import _EnvironmentPrints
+from ..tools import exponential_backoff
 
 try:
     import netCDF4
@@ -57,8 +58,7 @@ class Environment:
     Environment.datum : string
         The desired reference ellipsoid model, the following options are
         available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default
-        is "SIRGAS2000", then this model will be used if the user make some
-        typing mistake
+        is "SIRGAS2000".
     Environment.initial_east : float
         Launch site East UTM coordinate
     Environment.initial_north :  float
@@ -74,7 +74,7 @@ class Environment:
         Launch site E/W hemisphere
     Environment.elevation : float
         Launch site elevation.
-    Environment.date : datetime
+    Environment.datetime_date : datetime
         Date time of launch in UTC.
     Environment.local_date : datetime
         Date time of launch in the local time zone, defined by
@@ -276,49 +276,70 @@ def __init__(
         timezone="UTC",
         max_expected_height=80000.0,
     ):
-        """Initialize Environment class, saving launch rail length,
-        launch date, location coordinates and elevation. Note that
-        by default the standard atmosphere is loaded until another
+        """Initializes the Environment class, capturing essential parameters of
+        the launch site, including the launch date, geographical coordinates,
+        and elevation. This class is designed to calculate crucial variables
+        for the Flight simulation, such as atmospheric air pressure, density,
+        and gravitational acceleration.
+
+        Note that the default atmospheric model is the International Standard
+        Atmosphere as defined by ISO 2533 unless specified otherwise in
+        :meth:`Environment.set_atmospheric_model`.
 
         Parameters
         ----------
         gravity : int, float, callable, string, array, optional
             Surface gravitational acceleration. Positive values point the
-            acceleration down. If None, the Somigliana formula is used to
-        date : array, optional
-            Array of length 4, stating (year, month, day, hour (UTC))
-            of rocket launch. Must be given if a Forecast, Reanalysis
+            acceleration down. If None, the Somigliana formula is used.
+            See :meth:`Environment.set_gravity_model` for more information.
+        date : list or tuple, optional
+            List or tuple of length 4, stating (year, month, day, hour) in the
+            time zone of the parameter ``timezone``.
+            Alternatively, can be a ``datetime`` object specifying launch
+            date and time. The dates are stored as follows:
+
+            - :attr:`Environment.local_date`: Local time of launch in
+              the time zone specified by the parameter ``timezone``.
+
+            - :attr:`Environment.datetime_date`: UTC time of launch.
+
+            Must be given if a Forecast, Reanalysis
             or Ensemble, will be set as an atmospheric model.
+            Default is None.
+            See :meth:`Environment.set_date` for more information.
         latitude : float, optional
             Latitude in degrees (ranging from -90 to 90) of rocket
             launch location. Must be given if a Forecast, Reanalysis
             or Ensemble will be used as an atmospheric model or if
-            Open-Elevation will be used to compute elevation.
+            Open-Elevation will be used to compute elevation. Positive
+            values correspond to the North. Default value is 0, which
+            corresponds to the equator.
         longitude : float, optional
-            Longitude in degrees (ranging from -180 to 360) of rocket
+            Longitude in degrees (ranging from -180 to 180) of rocket
             launch location. Must be given if a Forecast, Reanalysis
             or Ensemble will be used as an atmospheric model or if
-            Open-Elevation will be used to compute elevation.
+            Open-Elevation will be used to compute elevation. Positive
+            values correspond to the East. Default value is 0, which
+            corresponds to the Greenwich Meridian.
         elevation : float, optional
             Elevation of launch site measured as height above sea
             level in meters. Alternatively, can be set as
             'Open-Elevation' which uses the Open-Elevation API to
             find elevation data. For this option, latitude and
             longitude must also be specified. Default value is 0.
-        datum : string
+        datum : string, optional
             The desired reference ellipsoidal model, the following options are
             available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default
-            is "SIRGAS2000", then this model will be used if the user make some
-            typing mistake.
+            is "SIRGAS2000".
         timezone : string, optional
             Name of the time zone. To see all time zones, import pytz and run
-            print(pytz.all_timezones). Default time zone is "UTC".
+            ``print(pytz.all_timezones)``. Default time zone is "UTC".
         max_expected_height : float, optional
             Maximum altitude in meters to keep weather data. The altitude must
             be above sea level (ASL). Especially useful for visualization.
             Can be altered as desired by doing `max_expected_height = number`.
             Depending on the atmospheric model, this value may be automatically
-            mofified.
+            modified.
 
         Returns
         -------
@@ -396,15 +417,57 @@ def set_date(self, date, timezone="UTC"):
 
         Parameters
         ----------
-        date : Datetime
-            Datetime object specifying launch date and time.
+        date : list, tuple, datetime
+            List or tuple of length 4, stating (year, month, day, hour) in the
+            time zone of the parameter ``timezone``. See Notes for more
+            information.
+            Alternatively, can be a ``datetime`` object specifying launch
+            date and time.
         timezone : string, optional
             Name of the time zone. To see all time zones, import pytz and run
-            print(pytz.all_timezones). Default time zone is "UTC".
+            ``print(pytz.all_timezones)``. Default time zone is "UTC".
 
         Returns
         -------
         None
+
+        Notes
+        -----
+        - If the ``date`` is given as a list or tuple, it should be in the same
+          time zone as specified by the ``timezone`` parameter. This local
+          time will be available in the attribute :attr:`Environment.local_date`
+          while the UTC time will be available in the attribute
+          :attr:`Environment.datetime_date`.
+
+        - If the ``date`` is given as a ``datetime`` object without a time zone,
+          it will be assumed to be in the same time zone as specified by the
+          ``timezone`` parameter. However, if the ``datetime`` object has a time
+          zone specified in its ``tzinfo`` attribute, the ``timezone``
+          parameter will be ignored.
+
+        Examples
+        --------
+
+        Let's set the launch date as an list:
+
+        >>> date = [2000, 1, 1, 13] # January 1st, 2000 at 13:00 UTC+1
+        >>> env = Environment()
+        >>> env.set_date(date, timezone="Europe/Rome")
+        >>> print(env.datetime_date) # Get UTC time
+        2000-01-01 12:00:00+00:00
+        >>> print(env.local_date)
+        2000-01-01 13:00:00+01:00
+
+        Now let's set the launch date as a ``datetime`` object:
+
+        >>> from datetime import datetime
+        >>> date = datetime(2000, 1, 1, 13, 0, 0)
+        >>> env = Environment()
+        >>> env.set_date(date, timezone="Europe/Rome")
+        >>> print(env.datetime_date) # Get UTC time
+        2000-01-01 12:00:00+00:00
+        >>> print(env.local_date)
+        2000-01-01 13:00:00+01:00
         """
         # Store date and configure time zone
         self.timezone = timezone
@@ -458,23 +521,66 @@ def set_location(self, latitude, longitude):
                 self.atmospheric_model_file, self.atmospheric_model_dict
             )
 
-        # Return None
-
-    def set_gravity_model(self, gravity):
-        """Sets the gravity model to be used in the simulation based on the
-        given user input to the gravity parameter.
+    def set_gravity_model(self, gravity=None):
+        """Defines the gravity model based on the given user input to the
+        gravity parameter. The gravity model is responsible for computing the
+        gravity acceleration at a given height above sea level in meters.
 
         Parameters
         ----------
-        gravity : None or Function source
-            If None, the Somigliana formula is used to compute the gravity
-            acceleration. Otherwise, the user can provide a Function object
-            representing the gravity model.
+        gravity : int, float, callable, string, list, optional
+            The gravitational acceleration in m/s² to be used in the
+            simulation, this value is positive when pointing downwards.
+            The input type can be one of the following:
+
+            - ``int`` or ``float``: The gravity acceleration is set as a\
+              constant function with respect to height;
+
+            - ``callable``: This callable should receive the height above\
+              sea level in meters and return the gravity acceleration;
+
+            - ``list``: The datapoints should be structured as\
+              ``[(h_i,g_i), ...]`` where ``h_i`` is the height above sea\
+              level in meters and ``g_i`` is the gravity acceleration in m/s²;
+
+            - ``string``: The string should correspond to a path to a CSV file\
+              containing the gravity acceleration data;
+
+            - ``None``: The Somigliana formula is used to compute the gravity\
+              acceleration.
+
+            This parameter is used as a :class:`Function` object source, check\
+            out the available input types for a more detailed explanation.
 
         Returns
         -------
         Function
             Function object representing the gravity model.
+
+        Notes
+        -----
+        This method **does not** set the gravity acceleration, it only returns
+        a :class:`Function` object representing the gravity model.
+
+        Examples
+        --------
+        Let's prepare a `Environment` object with a constant gravity
+        acceleration:
+
+        >>> g_0 = 9.80665
+        >>> env_cte_g = Environment(gravity=g_0)
+        >>> env_cte_g.gravity([0, 100, 1000])
+        [9.80665, 9.80665, 9.80665]
+
+        It's also possible to variate the gravity acceleration by defining
+        its function of height:
+
+        >>> R_t = 6371000
+        >>> g_func = lambda h : g_0 * (R_t / (R_t + h))**2
+        >>> env_var_g = Environment(gravity=g_func)
+        >>> g = env_var_g.gravity(1000)
+        >>> print(f"{g:.6f}")
+        9.803572
         """
         if gravity is None:
             return self.somigliana_gravity.set_discrete(
@@ -500,7 +606,7 @@ def max_expected_height(self, value):
 
     @funcify_method("height (m)", "gravity (m/s²)")
     def somigliana_gravity(self, height):
-        """Computes the gravity acceleration with the Somigliana formula.
+        """Computes the gravity acceleration with the Somigliana formula [1]_.
         An height correction is applied to the normal gravity that is
         accurate for heights used in aviation. The formula is based on the
         WGS84 ellipsoid, but is accurate for other reference ellipsoids.
@@ -514,6 +620,10 @@ def somigliana_gravity(self, height):
         -------
         Function
             Function object representing the gravity model.
+
+        References
+        ----------
+        .. [1] https://en.wikipedia.org/wiki/Theoretical_gravity#Somigliana_equation
         """
         a = 6378137.0  # semi_major_axis
         f = 1 / 298.257223563  # flattening_factor
@@ -571,18 +681,9 @@ def set_elevation(self, elevation="Open-Elevation"):
 
         #     self.elevation = elev
 
-        elif self.latitude != None and self.longitude != None:
-            try:
-                print("Fetching elevation from open-elevation.com...")
-                request_url = "https://api.open-elevation.com/api/v1/lookup?locations={:f},{:f}".format(
-                    self.latitude, self.longitude
-                )
-                response = requests.get(request_url)
-                results = response.json()["results"]
-                self.elevation = results[0]["elevation"]
-                print("Elevation received:", self.elevation)
-            except:
-                raise RuntimeError("Unable to reach Open-Elevation API servers.")
+        elif self.latitude is not None and self.longitude is not None:
+            self.elevation = self.__fetch_open_elevation()
+            print("Elevation received: ", self.elevation)
         else:
             raise ValueError(
                 "Latitude and longitude must be set to use"
@@ -1194,26 +1295,8 @@ def set_atmospheric_model(
                     "v_wind": "vgrdprs",
                 }
                 # Attempt to get latest forecast
-                time_attempt = datetime.utcnow()
-                success = False
-                attempt_count = 0
-                while not success and attempt_count < 10:
-                    time_attempt -= timedelta(hours=6 * attempt_count)
-                    file = "https://nomads.ncep.noaa.gov/dods/gens_bc/gens{:04d}{:02d}{:02d}/gep_all_{:02d}z".format(
-                        time_attempt.year,
-                        time_attempt.month,
-                        time_attempt.day,
-                        6 * (time_attempt.hour // 6),
-                    )
-                    try:
-                        self.process_ensemble(file, dictionary)
-                        success = True
-                    except OSError:
-                        attempt_count += 1
-                if not success:
-                    raise RuntimeError(
-                        "Unable to load latest weather data for GEFS through " + file
-                    )
+                self.__fetch_gefs_ensemble(dictionary)
+
             elif file == "CMC":
                 # Define dictionary
                 dictionary = {
@@ -1229,27 +1312,7 @@ def set_atmospheric_model(
                     "u_wind": "ugrdprs",
                     "v_wind": "vgrdprs",
                 }
-                # Attempt to get latest forecast
-                time_attempt = datetime.utcnow()
-                success = False
-                attempt_count = 0
-                while not success and attempt_count < 10:
-                    time_attempt -= timedelta(hours=12 * attempt_count)
-                    file = "https://nomads.ncep.noaa.gov/dods/cmcens/cmcens{:04d}{:02d}{:02d}/cmcens_all_{:02d}z".format(
-                        time_attempt.year,
-                        time_attempt.month,
-                        time_attempt.day,
-                        12 * (time_attempt.hour // 12),
-                    )
-                    try:
-                        self.process_ensemble(file, dictionary)
-                        success = True
-                    except OSError:
-                        attempt_count += 1
-                if not success:
-                    raise RuntimeError(
-                        "Unable to load latest weather data for CMC through " + file
-                    )
+                self.__fetch_cmc_ensemble(dictionary)
             # Process other forecasts or reanalysis
             else:
                 # Check if default dictionary was requested
@@ -1541,20 +1604,7 @@ def process_windy_atmosphere(self, model="ECMWF"):
             model.
         """
 
-        # Process the model string
-        model = model.lower()
-        if model[-1] == "u":  # case iconEu
-            model = "".join([model[:4], model[4].upper(), model[4 + 1 :]])
-        # Load data from Windy.com: json file
-        url = f"https://node.windy.com/forecast/meteogram/{model}/{self.latitude}/{self.longitude}/?step=undefined"
-        try:
-            response = requests.get(url).json()
-        except:
-            if model == "iconEu":
-                raise ValueError(
-                    "Could not get a valid response for Icon-EU from Windy. Check if the latitude and longitude coordinates set are inside Europe.",
-                )
-            raise
+        response = self.__fetch_atmospheric_data_from_windy(model)
 
         # Determine time index from model
         time_array = np.array(response["data"]["hours"])
@@ -1715,18 +1765,7 @@ def process_wyoming_sounding(self, file):
         None
         """
         # Request Wyoming Sounding from file url
-        response = requests.get(file)
-        if response.status_code != 200:
-            raise ImportError("Unable to load " + file + ".")
-        if len(re.findall("Can't get .+ Observations at", response.text)):
-            raise ValueError(
-                re.findall("Can't get .+ Observations at .+", response.text)[0]
-                + " Check station number and date."
-            )
-        if response.text == "Invalid OUTPUT: specified\n":
-            raise ValueError(
-                "Invalid OUTPUT: specified. Make sure the output is Text: List."
-            )
+        response = self.__fetch_wyoming_sounding(file)
 
         # Process Wyoming Sounding by finding data table and station info
         response_split_text = re.split("(<.{0,1}PRE>)", response.text)
@@ -1852,9 +1891,7 @@ def process_noaaruc_sounding(self, file):
         None
         """
         # Request NOAA Ruc Sounding from file url
-        response = requests.get(file)
-        if response.status_code != 200 or len(response.text) < 10:
-            raise ImportError("Unable to load " + file + ".")
+        response = self.__fetch_noaaruc_sounding(file)
 
         # Split response into lines
         lines = response.text.split("\n")
@@ -3429,6 +3466,110 @@ def set_earth_geometry(self, datum):
                 f"The reference system {datum} for Earth geometry " "is not recognized."
             )
 
+    # Auxiliary functions - Fetching Data from 3rd party APIs
+
+    @exponential_backoff(max_attempts=3, base_delay=1, max_delay=60)
+    def __fetch_open_elevation(self):
+        print("Fetching elevation from open-elevation.com...")
+        request_url = (
+            "https://api.open-elevation.com/api/v1/lookup?locations"
+            f"={self.latitude},{self.longitude}"
+        )
+        try:
+            response = requests.get(request_url)
+        except Exception as e:
+            raise RuntimeError("Unable to reach Open-Elevation API servers.")
+        results = response.json()["results"]
+        return results[0]["elevation"]
+
+    @exponential_backoff(max_attempts=5, base_delay=2, max_delay=60)
+    def __fetch_atmospheric_data_from_windy(self, model):
+        model = model.lower()
+        if model[-1] == "u":  # case iconEu
+            model = "".join([model[:4], model[4].upper(), model[4 + 1 :]])
+        url = (
+            f"https://node.windy.com/forecast/meteogram/{model}/"
+            f"{self.latitude}/{self.longitude}/?step=undefined"
+        )
+        try:
+            response = requests.get(url).json()
+        except Exception as e:
+            if model == "iconEu":
+                raise ValueError(
+                    "Could not get a valid response for Icon-EU from Windy. "
+                    "Check if the coordinates are set inside Europe."
+                )
+        return response
+
+    @exponential_backoff(max_attempts=5, base_delay=2, max_delay=60)
+    def __fetch_wyoming_sounding(self, file):
+        response = requests.get(file)
+        if response.status_code != 200:
+            raise ImportError(f"Unable to load {file}.")
+        if len(re.findall("Can't get .+ Observations at", response.text)):
+            raise ValueError(
+                re.findall("Can't get .+ Observations at .+", response.text)[0]
+                + " Check station number and date."
+            )
+        if response.text == "Invalid OUTPUT: specified\n":
+            raise ValueError(
+                "Invalid OUTPUT: specified. Make sure the output is Text: List."
+            )
+        return response
+
+    @exponential_backoff(max_attempts=5, base_delay=2, max_delay=60)
+    def __fetch_noaaruc_sounding(self, file):
+        response = requests.get(file)
+        if response.status_code != 200 or len(response.text) < 10:
+            raise ImportError("Unable to load " + file + ".")
+
+    @exponential_backoff(max_attempts=5, base_delay=2, max_delay=60)
+    def __fetch_gefs_ensemble(self, dictionary):
+        time_attempt = datetime.now(tz=timezone.utc)
+        success = False
+        attempt_count = 0
+        while not success and attempt_count < 10:
+            time_attempt -= timedelta(hours=6 * attempt_count)
+            file = (
+                f"https://nomads.ncep.noaa.gov/dods/gens_bc/gens"
+                f"{time_attempt.year:04d}{time_attempt.month:02d}"
+                f"{time_attempt.day:02d}/"
+                f"gep_all_{6 * (time_attempt.hour // 6):02d}z"
+            )
+            try:
+                self.process_ensemble(file, dictionary)
+                success = True
+            except OSError:
+                attempt_count += 1
+        if not success:
+            raise RuntimeError(
+                "Unable to load latest weather data for GEFS through " + file
+            )
+
+    @exponential_backoff(max_attempts=5, base_delay=2, max_delay=60)
+    def __fetch_cmc_ensemble(self, dictionary):
+        # Attempt to get latest forecast
+        time_attempt = datetime.now(tz=timezone.utc)
+        success = False
+        attempt_count = 0
+        while not success and attempt_count < 10:
+            time_attempt -= timedelta(hours=12 * attempt_count)
+            file = (
+                f"https://nomads.ncep.noaa.gov/dods/cmcens/"
+                f"cmcens{time_attempt.year:04d}{time_attempt.month:02d}"
+                f"{time_attempt.day:02d}/"
+                f"cmcens_all_{12 * (time_attempt.hour // 12):02d}z"
+            )
+            try:
+                self.process_ensemble(file, dictionary)
+                success = True
+            except OSError:
+                attempt_count += 1
+        if not success:
+            raise RuntimeError(
+                "Unable to load latest weather data for CMC through " + file
+            )
+
     # Auxiliary functions - Geodesic Coordinates
 
     @staticmethod
diff --git a/rocketpy/environment/environment_analysis.py b/rocketpy/environment/environment_analysis.py
index c15b32551..da6fde364 100644
--- a/rocketpy/environment/environment_analysis.py
+++ b/rocketpy/environment/environment_analysis.py
@@ -3,6 +3,7 @@
 import datetime
 import json
 from collections import defaultdict
+from functools import cached_property
 
 import netCDF4
 import numpy as np
@@ -22,11 +23,6 @@
 from ..units import convert_units
 from .environment import Environment
 
-try:
-    from functools import cached_property
-except ImportError:
-    from ..tools import cached_property
-
 # TODO: the average_wind_speed_profile_by_hour and similar methods could be more abstract than currently are
 
 
@@ -441,7 +437,7 @@ def __localize_input_dates(self):
 
     def __find_preferred_timezone(self):
         if self.preferred_timezone is None:
-            # Use local timezone based on lat lon pair
+            # Use local time zone based on lat lon pair
             try:
                 timezonefinder = import_optional_dependency("timezonefinder")
                 tf = timezonefinder.TimezoneFinder()
diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py
index cefed044d..eda903bcc 100644
--- a/rocketpy/mathutils/function.py
+++ b/rocketpy/mathutils/function.py
@@ -7,6 +7,7 @@
 import warnings
 from collections.abc import Iterable
 from copy import deepcopy
+from functools import cached_property
 from inspect import signature
 from pathlib import Path
 
@@ -14,11 +15,6 @@
 import numpy as np
 from scipy import integrate, linalg, optimize
 
-try:
-    from functools import cached_property
-except ImportError:
-    from ..tools import cached_property
-
 NUMERICAL_TYPES = (float, int, complex, np.ndarray, np.integer, np.floating)
 
 
diff --git a/rocketpy/mathutils/vector_matrix.py b/rocketpy/mathutils/vector_matrix.py
index 7a658de04..9c3efe616 100644
--- a/rocketpy/mathutils/vector_matrix.py
+++ b/rocketpy/mathutils/vector_matrix.py
@@ -1,7 +1,8 @@
 from cmath import isclose
+from functools import cached_property
 from itertools import product
 
-from rocketpy.tools import cached_property, euler_to_quaternions
+from rocketpy.tools import euler_to_quaternions
 
 
 class Vector:
diff --git a/rocketpy/motors/hybrid_motor.py b/rocketpy/motors/hybrid_motor.py
index 6f0849cd0..557333fe7 100644
--- a/rocketpy/motors/hybrid_motor.py
+++ b/rocketpy/motors/hybrid_motor.py
@@ -1,3 +1,5 @@
+from functools import cached_property
+
 from rocketpy.tools import parallel_axis_theorem_from_com
 
 from ..mathutils.function import Function, funcify_method, reset_funcified_methods
@@ -7,11 +9,6 @@
 from .motor import Motor
 from .solid_motor import SolidMotor
 
-try:
-    from functools import cached_property
-except ImportError:
-    from ..tools import cached_property
-
 
 class HybridMotor(Motor):
     """Class to specify characteristics and useful operations for Hybrid
diff --git a/rocketpy/motors/liquid_motor.py b/rocketpy/motors/liquid_motor.py
index 01f728473..7314e11ba 100644
--- a/rocketpy/motors/liquid_motor.py
+++ b/rocketpy/motors/liquid_motor.py
@@ -1,4 +1,4 @@
-import warnings
+from functools import cached_property
 
 import numpy as np
 
@@ -13,11 +13,6 @@
 from ..prints.liquid_motor_prints import _LiquidMotorPrints
 from .motor import Motor
 
-try:
-    from functools import cached_property
-except ImportError:
-    from ..tools import cached_property
-
 
 class LiquidMotor(Motor):
     """Class to specify characteristics and useful operations for Liquid
diff --git a/rocketpy/motors/motor.py b/rocketpy/motors/motor.py
index 3834f4a15..9429da88e 100644
--- a/rocketpy/motors/motor.py
+++ b/rocketpy/motors/motor.py
@@ -1,6 +1,7 @@
 import re
 import warnings
 from abc import ABC, abstractmethod
+from functools import cached_property
 
 import numpy as np
 
@@ -9,11 +10,6 @@
 from ..prints.motor_prints import _MotorPrints
 from ..tools import parallel_axis_theorem_from_com, tuple_handler
 
-try:
-    from functools import cached_property
-except ImportError:
-    from ..tools import cached_property
-
 
 class Motor(ABC):
     """Abstract class to specify characteristics and useful operations for
diff --git a/rocketpy/motors/solid_motor.py b/rocketpy/motors/solid_motor.py
index 8b1c2362e..db3527a95 100644
--- a/rocketpy/motors/solid_motor.py
+++ b/rocketpy/motors/solid_motor.py
@@ -1,3 +1,5 @@
+from functools import cached_property
+
 import numpy as np
 from scipy import integrate
 
@@ -6,11 +8,6 @@
 from ..prints.solid_motor_prints import _SolidMotorPrints
 from .motor import Motor
 
-try:
-    from functools import cached_property
-except ImportError:
-    from ..tools import cached_property
-
 
 class SolidMotor(Motor):
     """Class to specify characteristics and useful operations for solid motors.
diff --git a/rocketpy/motors/tank_geometry.py b/rocketpy/motors/tank_geometry.py
index f1940cbea..2eb7bd27e 100644
--- a/rocketpy/motors/tank_geometry.py
+++ b/rocketpy/motors/tank_geometry.py
@@ -11,10 +11,7 @@
 
     cache = lru_cache(maxsize=None)
 
-try:
-    from functools import cached_property
-except ImportError:
-    from ..tools import cached_property
+from functools import cached_property
 
 
 class TankGeometry:
diff --git a/rocketpy/plots/flight_plots.py b/rocketpy/plots/flight_plots.py
index 4ae5141c9..21266a1f3 100644
--- a/rocketpy/plots/flight_plots.py
+++ b/rocketpy/plots/flight_plots.py
@@ -1,11 +1,8 @@
+from functools import cached_property
+
 import matplotlib.pyplot as plt
 import numpy as np
 
-try:
-    from functools import cached_property
-except ImportError:
-    from ..tools import cached_property
-
 
 class _FlightPlots:
     """Class that holds plot methods for Flight class.
diff --git a/rocketpy/plots/rocket_plots.py b/rocketpy/plots/rocket_plots.py
index e801da3c7..0d7b5b130 100644
--- a/rocketpy/plots/rocket_plots.py
+++ b/rocketpy/plots/rocket_plots.py
@@ -230,6 +230,8 @@ def draw(self, vis_args=None, plane="xz"):
         plt.show()
 
     def _draw_aerodynamic_surfaces(self, ax, vis_args):
+        """Draws the aerodynamic surfaces and saves the position of the points
+        of interest for the tubes."""
         # List of drawn surfaces with the position of points of interest
         # and the radius of the rocket at that point
         drawn_surfaces = []
@@ -250,6 +252,8 @@ def _draw_aerodynamic_surfaces(self, ax, vis_args):
         return drawn_surfaces
 
     def _draw_nose_cone(self, ax, surface, position, drawn_surfaces, vis_args):
+        """Draws the nosecone and saves the position of the points of interest
+        for the tubes."""
         x_nosecone = -self.rocket._csys * surface.shape_vec[0] + position
         y_nosecone = surface.shape_vec[1]
         ax.plot(
@@ -277,6 +281,8 @@ def _draw_nose_cone(self, ax, surface, position, drawn_surfaces, vis_args):
         )
 
     def _draw_tail(self, ax, surface, position, drawn_surfaces, vis_args):
+        """Draws the tail and saves the position of the points of interest
+        for the tubes."""
         x_tail = -self.rocket._csys * surface.shape_vec[0] + position
         y_tail = surface.shape_vec[1]
         ax.plot(
@@ -302,13 +308,30 @@ def _draw_tail(self, ax, surface, position, drawn_surfaces, vis_args):
         drawn_surfaces.append((surface, position, surface.bottom_radius, x_tail[-1]))
 
     def _draw_fins(self, ax, surface, position, drawn_surfaces, vis_args):
+        """Draws the fins and saves the position of the points of interest
+        for the tubes."""
         num_fins = surface.n
         x_fin = -self.rocket._csys * surface.shape_vec[0] + position
         y_fin = surface.shape_vec[1] + surface.rocket_radius
         rotation_angles = [2 * np.pi * i / num_fins for i in range(num_fins)]
 
         for angle in rotation_angles:
-            x_rotated, y_rotated = self._rotate_points(x_fin, y_fin, angle)
+            # Create a rotation matrix for the current angle around the x-axis
+            rotation_matrix = np.array([[1, 0], [0, np.cos(angle)]])
+
+            # Apply the rotation to the original fin points
+            rotated_points_2d = np.dot(rotation_matrix, np.vstack((x_fin, y_fin)))
+
+            # Extract x and y coordinates of the rotated points
+            x_rotated, y_rotated = rotated_points_2d
+
+            # Project points above the XY plane back into the XY plane (set z-coordinate to 0)
+            x_rotated = np.where(
+                rotated_points_2d[1] > 0, rotated_points_2d[0], x_rotated
+            )
+            y_rotated = np.where(
+                rotated_points_2d[1] > 0, rotated_points_2d[1], y_rotated
+            )
             ax.plot(
                 x_rotated,
                 y_rotated,
@@ -318,22 +341,8 @@ def _draw_fins(self, ax, surface, position, drawn_surfaces, vis_args):
 
         drawn_surfaces.append((surface, position, surface.rocket_radius, x_rotated[-1]))
 
-    def _rotate_points(self, x_fin, y_fin, angle):
-        # Create a rotation matrix for the current angle around the x-axis
-        rotation_matrix = np.array([[1, 0], [0, np.cos(angle)]])
-
-        # Apply the rotation to the original fin points
-        rotated_points_2d = np.dot(rotation_matrix, np.vstack((x_fin, y_fin)))
-
-        # Extract x and y coordinates of the rotated points
-        x_rotated, y_rotated = rotated_points_2d
-
-        # Project points above the XY plane back into the XY plane (set z-coordinate to 0)
-        x_rotated = np.where(rotated_points_2d[1] > 0, rotated_points_2d[0], x_rotated)
-        y_rotated = np.where(rotated_points_2d[1] > 0, rotated_points_2d[1], y_rotated)
-        return x_rotated, y_rotated
-
     def _draw_tubes(self, ax, drawn_surfaces, vis_args):
+        """Draws the tubes between the aerodynamic surfaces."""
         for i, d_surface in enumerate(drawn_surfaces):
             # Draw the tubes, from the end of the first surface to the beginning
             # of the next surface, with the radius of the rocket at that point
@@ -373,6 +382,7 @@ def _draw_tubes(self, ax, drawn_surfaces, vis_args):
         return radius, last_x
 
     def _draw_motor(self, last_radius, last_x, ax, vis_args):
+        """Draws the motor from motor patches"""
         total_csys = self.rocket._csys * self.rocket.motor._csys
         nozzle_position = (
             self.rocket.motor_position + self.rocket.motor.nozzle_position * total_csys
@@ -400,6 +410,7 @@ def _draw_motor(self, last_radius, last_x, ax, vis_args):
         self._draw_nozzle_tube(last_radius, last_x, nozzle_position, ax, vis_args)
 
     def _generate_motor_patches(self, total_csys, ax, vis_args):
+        """Generates motor patches for drawing"""
         motor_patches = []
 
         if isinstance(self.rocket.motor, SolidMotor):
@@ -478,6 +489,7 @@ def _generate_motor_patches(self, total_csys, ax, vis_args):
         return motor_patches
 
     def _draw_nozzle_tube(self, last_radius, last_x, nozzle_position, ax, vis_args):
+        """Draws the tube from the last surface to the nozzle position."""
         # Check if nozzle is beyond the last surface, if so draw a tube
         # to it, with the radius of the last surface
         if self.rocket._csys == 1:
@@ -518,6 +530,7 @@ def _draw_nozzle_tube(self, last_radius, last_x, nozzle_position, ax, vis_args):
                 )
 
     def _draw_rail_buttons(self, ax, vis_args):
+        """Draws the rail buttons of the rocket."""
         try:
             buttons, pos = self.rocket.rail_buttons[0]
             lower = pos
@@ -532,6 +545,7 @@ def _draw_rail_buttons(self, ax, vis_args):
             pass
 
     def _draw_center_of_mass_and_pressure(self, ax):
+        """Draws the center of mass and center of pressure of the rocket."""
         # Draw center of mass and center of pressure
         cm = self.rocket.center_of_mass(0)
         ax.scatter(cm, 0, color="#1565c0", label="Center of Mass", s=10)
diff --git a/rocketpy/rocket/aero_surface.py b/rocketpy/rocket/aero_surface.py
index c5d154f3e..9649d1bf0 100644
--- a/rocketpy/rocket/aero_surface.py
+++ b/rocketpy/rocket/aero_surface.py
@@ -1,6 +1,5 @@
 import warnings
 from abc import ABC, abstractmethod
-from functools import cached_property
 
 import numpy as np
 from scipy.optimize import fsolve
@@ -1978,8 +1977,9 @@ def __init__(
             brakes drag coefficient will be used for the entire rocket instead.
             Default is False.
         deployment_level : float, optional
-            Current deployment level, ranging from 0 to 1. Deployment level is the
-            fraction of the total airbrake area that is Deployment. Default is 0.
+            Initial deployment level, ranging from 0 to 1. Deployment level is
+            the fraction of the total airbrake area that is Deployment. Default
+            is 0.
         name : str, optional
             Name of the air brakes. Default is "AirBrakes".
 
@@ -1997,6 +1997,7 @@ def __init__(
         self.reference_area = reference_area
         self.clamp = clamp
         self.override_rocket_drag = override_rocket_drag
+        self.initial_deployment_level = deployment_level
         self.deployment_level = deployment_level
         self.prints = _AirBrakesPrints(self)
         self.plots = _AirBrakesPlots(self)
@@ -2023,6 +2024,12 @@ def deployment_level(self, value):
                 )
         self._deployment_level = value
 
+    def _reset(self):
+        """Resets the air brakes to their initial state. This is ran at the
+        beginning of each simulation to ensure the air brakes are in the correct
+        state."""
+        self.deployment_level = self.initial_deployment_level
+
     def evaluate_center_of_pressure(self):
         """Evaluates the center of pressure of the aerodynamic surface in local
         coordinates.
diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py
index 9c179557e..45fce0cd5 100644
--- a/rocketpy/rocket/rocket.py
+++ b/rocketpy/rocket/rocket.py
@@ -82,6 +82,9 @@ class Rocket:
         Function of time expressing the total mass of the rocket,
         defined as the sum of the propellant mass and the rocket
         mass without propellant.
+    Rocket.total_mass_flow_rate : Function
+        Time derivative of rocket's total mass in kg/s as a function
+        of time as obtained by the thrust source of the added motor.
     Rocket.thrust_to_weight : Function
         Function of time expressing the motor thrust force divided by rocket
         weight. The gravitational acceleration is assumed as 9.80665 m/s^2.
@@ -760,6 +763,7 @@ def add_motor(self, motor, position):
             self.motor.center_of_dry_mass_position * _ + self.motor_position
         )
         self.nozzle_position = self.motor.nozzle_position * _ + self.motor_position
+        self.total_mass_flow_rate = self.motor.total_mass_flow_rate
         self.evaluate_dry_mass()
         self.evaluate_total_mass()
         self.evaluate_center_of_dry_mass()
diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index f82b63b1c..8aa3f4cd6 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -611,6 +611,9 @@ def __init__(
         self.name = name
         self.equations_of_motion = equations_of_motion
 
+        # Controller initialization
+        self.__init_controllers()
+
         # Flight initialization
         self.__init_post_process_variables()
         self.__init_solution_monitors()
@@ -1107,8 +1110,14 @@ def __init__(
                                             [self.t, parachute]
                                         )
 
+                    # If controlled flight, post process must be done on sim time
+                    if self._controllers:
+                        phase.derivative(self.t, self.y_sol, post_processing=True)
+
         self.t_final = self.t
         self._calculate_pressure_signal()
+        if self._controllers:
+            self.__cache_post_process_variables()
         if verbose:
             print("Simulation Completed at Time: {:3.4f} s".format(self.t))
 
@@ -1120,6 +1129,25 @@ def __init_post_process_variables(self):
         self._bearing = Function(0)
         self._latitude = Function(0)
         self._longitude = Function(0)
+        # Initialize state derivatives, force and atmospheric arrays
+        self.ax_list = []
+        self.ay_list = []
+        self.az_list = []
+        self.alpha1_list = []
+        self.alpha2_list = []
+        self.alpha3_list = []
+        self.R1_list = []
+        self.R2_list = []
+        self.R3_list = []
+        self.M1_list = []
+        self.M2_list = []
+        self.M3_list = []
+        self.pressure_list = []
+        self.density_list = []
+        self.dynamic_viscosity_list = []
+        self.speed_of_sound_list = []
+        self.wind_velocity_x_list = []
+        self.wind_velocity_y_list = []
 
     def __init_solution_monitors(self):
         # Initialize solution monitors
@@ -1192,6 +1220,11 @@ def __init_flight_state(self):
             self.out_of_rail_time = self.initial_solution[0]
             self.out_of_rail_time_index = 0
             self.initial_derivative = self.u_dot_generalized
+        if self._controllers:
+            # Handle post process during simulation, get initial accel/forces
+            self.initial_derivative(
+                self.t_initial, self.initial_solution[1:], post_processing=True
+            )
 
     def __init_solver_monitors(self):
         # Initialize solver monitors
@@ -1212,10 +1245,41 @@ def __init_equations_of_motion(self):
         if self.equations_of_motion == "solid_propulsion":
             self.u_dot_generalized = self.u_dot
 
-    def __init_equations_of_motion(self):
-        """Initialize equations of motion."""
-        if self.equations_of_motion == "solid_propulsion":
-            self.u_dot_generalized = self.u_dot
+    def __init_controllers(self):
+        """Initialize controllers"""
+        self._controllers = self.rocket._controllers[:]
+        if self._controllers:
+            if self.time_overshoot == True:
+                self.time_overshoot = False
+                warnings.warn(
+                    "time_overshoot has been set to False due to the presence of controllers. "
+                )
+            # reset controllable object to initial state (only airbrakes for now)
+            for air_brakes in self.rocket.air_brakes:
+                air_brakes._reset()
+
+    def __cache_post_process_variables(self):
+        """Cache post-process variables for simulations with controllers."""
+        self.__retrieve_arrays = [
+            self.ax_list,
+            self.ay_list,
+            self.az_list,
+            self.alpha1_list,
+            self.alpha2_list,
+            self.alpha3_list,
+            self.R1_list,
+            self.R2_list,
+            self.R3_list,
+            self.M1_list,
+            self.M2_list,
+            self.M3_list,
+            self.pressure_list,
+            self.density_list,
+            self.dynamic_viscosity_list,
+            self.speed_of_sound_list,
+            self.wind_velocity_x_list,
+            self.wind_velocity_y_list,
+        ]
 
     @cached_property
     def effective_1rl(self):
@@ -1292,11 +1356,6 @@ def udot_rail1(self, t, u, post_processing=False):
             e0dot, e1dot, e2dot, e3dot, alpha1, alpha2, alpha3].
 
         """
-        # Check if post processing mode is on
-        if post_processing:
-            # Use u_dot post processing code
-            return self.u_dot_generalized(t, u, True)
-
         # Retrieve integration data
         x, y, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u
 
@@ -1327,6 +1386,17 @@ def udot_rail1(self, t, u, post_processing=False):
         else:
             ax, ay, az = 0, 0, 0
 
+        if post_processing:
+            # Use u_dot post processing code for forces, moments and env data
+            self.u_dot_generalized(t, u, post_processing=True)
+            # Save feasible accelerations
+            self.ax_list[-1] = [t, ax]
+            self.ay_list[-1] = [t, ay]
+            self.az_list[-1] = [t, az]
+            self.alpha1_list[-1] = [t, 0]
+            self.alpha2_list[-1] = [t, 0]
+            self.alpha3_list[-1] = [t, 0]
+
         return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0]
 
     def udot_rail2(self, t, u, post_processing=False):
@@ -1616,6 +1686,13 @@ def u_dot(self, t, u, post_processing=False):
         ]
 
         if post_processing:
+            # Accelerations
+            self.ax_list.append([t, ax])
+            self.ay_list.append([t, ay])
+            self.az_list.append([t, az])
+            self.alpha1_list.append([t, alpha1])
+            self.alpha2_list.append([t, alpha2])
+            self.alpha3_list.append([t, alpha3])
             # Dynamics variables
             self.R1_list.append([t, R1])
             self.R2_list.append([t, R2])
@@ -1676,8 +1753,8 @@ def u_dot_generalized(self, t, u, post_processing=False):
         # Retrieve necessary quantities
         rho = self.env.density.get_value_opt(z)
         total_mass = self.rocket.total_mass.get_value_opt(t)
-        total_mass_dot = self.rocket.total_mass.differentiate(t)
-        total_mass_ddot = self.rocket.total_mass.differentiate(t, order=2)
+        total_mass_dot = self.rocket.total_mass_flow_rate.get_value_opt(t)
+        total_mass_ddot = self.rocket.total_mass_flow_rate.differentiate(t)
         ## CM position vector and time derivatives relative to CDM in body frame
         r_CM_z = (
             -1
@@ -1891,6 +1968,13 @@ def u_dot_generalized(self, t, u, post_processing=False):
         u_dot = [*r_dot, *v_dot, *e_dot, *w_dot]
 
         if post_processing:
+            # Accelerations
+            self.ax_list.append([t, v_dot[0]])
+            self.ay_list.append([t, v_dot[1]])
+            self.az_list.append([t, v_dot[2]])
+            self.alpha1_list.append([t, w_dot[0]])
+            self.alpha2_list.append([t, w_dot[1]])
+            self.alpha3_list.append([t, w_dot[2]])
             # Dynamics variables
             self.R1_list.append([t, R1])
             self.R2_list.append([t, R2])
@@ -1977,6 +2061,13 @@ def u_dot_parachute(self, t, u, post_processing=False):
         u_dot = [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0]
 
         if post_processing:
+            # Accelerations
+            self.ax_list.append([t, ax])
+            self.ay_list.append([t, ay])
+            self.az_list.append([t, az])
+            self.alpha1_list.append([t, 0])
+            self.alpha2_list.append([t, 0])
+            self.alpha3_list.append([t, 0])
             # Dynamics variables
             self.R1_list.append([t, Dx])
             self.R2_list.append([t, Dy])
@@ -1985,13 +2076,20 @@ def u_dot_parachute(self, t, u, post_processing=False):
             self.M2_list.append([t, 0])
             self.M3_list.append([t, 0])
             # Atmospheric Conditions
-            self.wind_velocity_x_list.append([t, self.env.wind_velocity_x(z)])
-            self.wind_velocity_y_list.append([t, self.env.wind_velocity_y(z)])
-            self.density_list.append([t, self.env.density(z)])
-            self.dynamic_viscosity_list.append([t, self.env.dynamic_viscosity(z)])
-            self.pressure_list.append([t, self.env.pressure(z)])
-            self.speed_of_sound_list.append([t, self.env.speed_of_sound(z)])
-
+            self.wind_velocity_x_list.append(
+                [t, self.env.wind_velocity_x.get_value_opt(z)]
+            )
+            self.wind_velocity_y_list.append(
+                [t, self.env.wind_velocity_y.get_value_opt(z)]
+            )
+            self.density_list.append([t, self.env.density.get_value_opt(z)])
+            self.dynamic_viscosity_list.append(
+                [t, self.env.dynamic_viscosity.get_value_opt(z)]
+            )
+            self.pressure_list.append([t, self.env.pressure.get_value_opt(z)])
+            self.speed_of_sound_list.append(
+                [t, self.env.speed_of_sound.get_value_opt(z)]
+            )
         return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0]
 
     @cached_property
@@ -2824,31 +2922,20 @@ def longitude(self):
         return np.column_stack((self.time, longitude))
 
     @cached_property
-    def retrieve_acceleration_arrays(self):
-        """Retrieve acceleration arrays from the integration scheme
-
-        Parameters
-        ----------
+    def __retrieve_arrays(self):
+        """post processing function to retrieve arrays from the integration
+        scheme and store them in lists for further analysis.
 
         Returns
         -------
-        ax: list
-            acceleration in x direction
-        ay: list
-            acceleration in y direction
-        az: list
-            acceleration in z direction
-        alpha1: list
-            angular acceleration in x direction
-        alpha2: list
-            angular acceleration in y direction
-        alpha3: list
-            angular acceleration in z direction
+        temp_values: list
+            List containing the following arrays: ``ax`` , ``ay`` , ``az`` ,
+            ``alpha1`` , ``alpha2`` , ``alpha3`` , ``R1`` , ``R2`` , ``R3`` ,
+            ``M1`` , ``M2`` , ``M3`` , ``pressure`` , ``density`` ,
+            ``dynamic_viscosity`` , ``speed_of_sound`` , ``wind_velocity_x`` ,
+            ``wind_velocity_y``.
         """
-        # Initialize acceleration arrays
-        ax, ay, az = [[0, 0]], [[0, 0]], [[0, 0]]
-        alpha1, alpha2, alpha3 = [[0, 0]], [[0, 0]], [[0, 0]]
-        # Go through each time step and calculate accelerations
+        # Go through each time step and calculate forces and atmospheric values
         # Get flight phases
         for phase_index, phase in self.time_iterator(self.FlightPhases):
             init_time = phase.t
@@ -2857,23 +2944,60 @@ def retrieve_acceleration_arrays(self):
             # Call callback functions
             for callback in phase.callbacks:
                 callback(self)
-            # Loop through time steps in flight phase
-            for step in self.solution:  # Can be optimized
-                if init_time < step[0] <= final_time:
-                    # Get derivatives
-                    u_dot = current_derivative(step[0], step[1:])
-                    # Get accelerations
-                    ax_value, ay_value, az_value = u_dot[3:6]
-                    alpha1_value, alpha2_value, alpha3_value = u_dot[10:]
-                    # Save accelerations
-                    ax.append([step[0], ax_value])
-                    ay.append([step[0], ay_value])
-                    az.append([step[0], az_value])
-                    alpha1.append([step[0], alpha1_value])
-                    alpha2.append([step[0], alpha2_value])
-                    alpha3.append([step[0], alpha3_value])
-
-        return ax, ay, az, alpha1, alpha2, alpha3
+            # find index of initial and final time of phase in solution array
+            init_time_index = find_closest(self.time, init_time)
+            final_time_index = find_closest(self.time, final_time) + 1
+            # Loop through time steps solution array
+            for step in self.solution[init_time_index:final_time_index]:
+                if init_time != step[0] or (
+                    init_time == self.t_initial and step[0] == self.t_initial
+                ):
+                    # Call derivatives in post processing mode
+                    current_derivative(step[0], step[1:], post_processing=True)
+
+        temp_values = [
+            self.ax_list,
+            self.ay_list,
+            self.az_list,
+            self.alpha1_list,
+            self.alpha2_list,
+            self.alpha3_list,
+            self.R1_list,
+            self.R2_list,
+            self.R3_list,
+            self.M1_list,
+            self.M2_list,
+            self.M3_list,
+            self.pressure_list,
+            self.density_list,
+            self.dynamic_viscosity_list,
+            self.speed_of_sound_list,
+            self.wind_velocity_x_list,
+            self.wind_velocity_y_list,
+        ]
+
+        return temp_values
+
+    @cached_property
+    def retrieve_acceleration_arrays(self):
+        """Retrieve acceleration arrays from the integration scheme
+
+        Returns
+        -------
+        ax_list: list
+            acceleration in x direction
+        ay_list: list
+            acceleration in y direction
+        az_list: list
+            acceleration in z direction
+        alpha1_list: list
+            angular acceleration in x direction
+        alpha2_list: list
+            angular acceleration in y direction
+        alpha3_list: list
+            angular acceleration in z direction
+        """
+        return self.__retrieve_arrays[:6]
 
     @cached_property
     def retrieve_temporary_values_arrays(self):
@@ -2910,54 +3034,7 @@ def retrieve_temporary_values_arrays(self):
         self.wind_velocity_y_list: list
             Wind velocity in y direction at each time step.
         """
-
-        # Initialize force and atmospheric arrays
-        self.R1_list = []
-        self.R2_list = []
-        self.R3_list = []
-        self.M1_list = []
-        self.M2_list = []
-        self.M3_list = []
-        self.pressure_list = []
-        self.density_list = []
-        self.dynamic_viscosity_list = []
-        self.speed_of_sound_list = []
-        self.wind_velocity_x_list = []
-        self.wind_velocity_y_list = []
-
-        # Go through each time step and calculate forces and atmospheric values
-        # Get flight phases
-        for phase_index, phase in self.time_iterator(self.FlightPhases):
-            init_time = phase.t
-            final_time = self.FlightPhases[phase_index + 1].t
-            current_derivative = phase.derivative
-            # Call callback functions
-            for callback in phase.callbacks:
-                callback(self)
-            # Loop through time steps in flight phase
-            for step in self.solution:  # Can be optimized
-                if init_time < step[0] <= final_time or (
-                    init_time == self.t_initial and step[0] == self.t_initial
-                ):
-                    # Call derivatives in post processing mode
-                    u_dot = current_derivative(step[0], step[1:], post_processing=True)
-
-        temporary_values = [
-            self.R1_list,
-            self.R2_list,
-            self.R3_list,
-            self.M1_list,
-            self.M2_list,
-            self.M3_list,
-            self.pressure_list,
-            self.density_list,
-            self.dynamic_viscosity_list,
-            self.speed_of_sound_list,
-            self.wind_velocity_x_list,
-            self.wind_velocity_y_list,
-        ]
-
-        return temporary_values
+        return self.__retrieve_arrays[6:]
 
     def get_controller_observed_variables(self):
         """Retrieve the observed variables related to air brakes from the
diff --git a/rocketpy/tools.py b/rocketpy/tools.py
index 5666fdbdf..edea8aada 100644
--- a/rocketpy/tools.py
+++ b/rocketpy/tools.py
@@ -1,6 +1,8 @@
+import functools
 import importlib
 import importlib.metadata
 import re
+import time
 from bisect import bisect_left
 
 import numpy as np
@@ -8,36 +10,10 @@
 from cftime import num2pydate
 from packaging import version as packaging_version
 
-_NOT_FOUND = object()
-
 # Mapping of module name and the name of the package that should be installed
 INSTALL_MAPPING = {"IPython": "ipython"}
 
 
-class cached_property:
-    def __init__(self, func):
-        self.func = func
-        self.attrname = None
-        self.__doc__ = func.__doc__
-
-    def __set_name__(self, owner, name):
-        self.attrname = name
-
-    def __get__(self, instance, owner=None):
-        if instance is None:
-            return self
-        if self.attrname is None:
-            raise TypeError(
-                "Cannot use cached_property instance without calling __set_name__ on it."
-            )
-        cache = instance.__dict__
-        val = cache.get(self.attrname, _NOT_FOUND)
-        if val is _NOT_FOUND:
-            val = self.func(instance)
-            cache[self.attrname] = val
-        return val
-
-
 def tuple_handler(value):
     """Transforms the input value into a tuple that
     represents a range. If the input is an int or float,
@@ -153,7 +129,7 @@ def time_num_to_date_string(time_num, units, timezone, calendar="gregorian"):
     """Convert time number (usually hours before a certain date) into two
     strings: one for the date (example: 2022.04.31) and one for the hour
     (example: 14). See cftime.num2date for details on units and calendar.
-    Automatically converts time number from UTC to local timezone based on
+    Automatically converts time number from UTC to local time zone based on
     lat, lon coordinates. This function was created originally for the
     EnvironmentAnalysis class.
 
@@ -382,6 +358,25 @@ def check_requirement_version(module_name, version):
     return True
 
 
+def exponential_backoff(max_attempts, base_delay=1, max_delay=60):
+    def decorator(func):
+        @functools.wraps(func)
+        def wrapper(*args, **kwargs):
+            delay = base_delay
+            for i in range(max_attempts):
+                try:
+                    return func(*args, **kwargs)
+                except Exception as e:
+                    if i == max_attempts - 1:
+                        raise e from None
+                    delay = min(delay * 2, max_delay)
+                    time.sleep(delay)
+
+        return wrapper
+
+    return decorator
+
+
 def parallel_axis_theorem_from_com(com_inertia_moment, mass, distance):
     """Calculates the moment of inertia of a object relative to a new axis using
     the parallel axis theorem. The new axis is parallel to and at a distance
diff --git a/tests/fixtures/environment/environment_fixtures.py b/tests/fixtures/environment/environment_fixtures.py
index 8949f9973..851be3203 100644
--- a/tests/fixtures/environment/environment_fixtures.py
+++ b/tests/fixtures/environment/environment_fixtures.py
@@ -1,6 +1,7 @@
 from datetime import datetime, timedelta
 
 import pytest
+
 from rocketpy import Environment, EnvironmentAnalysis
 
 
@@ -54,8 +55,8 @@ def env_analysis():
     EnvironmentAnalysis
     """
     env_analysis = EnvironmentAnalysis(
-        start_date=datetime.datetime(2019, 10, 23),
-        end_date=datetime.datetime(2021, 10, 23),
+        start_date=datetime(2019, 10, 23),
+        end_date=datetime(2021, 10, 23),
         latitude=39.3897,
         longitude=-8.28896388889,
         start_hour=6,
diff --git a/tests/test_environment.py b/tests/test_environment.py
index 5fa0e2c45..7349d512b 100644
--- a/tests/test_environment.py
+++ b/tests/test_environment.py
@@ -1,5 +1,4 @@
 import datetime
-import time
 from unittest.mock import patch
 
 import pytest
@@ -64,13 +63,8 @@ def test_wyoming_sounding_atmosphere(mock_show, example_plain_env):
     #        "file" option, instead of receiving the URL as a string.
     URL = "http://weather.uwyo.edu/cgi-bin/sounding?region=samer&TYPE=TEXT%3ALIST&YEAR=2019&MONTH=02&FROM=0500&TO=0512&STNM=83779"
     # give it at least 5 times to try to download the file
-    for i in range(5):
-        try:
-            example_plain_env.set_atmospheric_model(type="wyoming_sounding", file=URL)
-            break
-        except:
-            time.sleep(1)  # wait 1 second before trying again
-            pass
+    example_plain_env.set_atmospheric_model(type="wyoming_sounding", file=URL)
+
     assert example_plain_env.all_info() == None
     assert abs(example_plain_env.pressure(0) - 93600.0) < 1e-8
     assert (
diff --git a/tests/test_flight.py b/tests/test_flight.py
index db882e185..4fb4036eb 100644
--- a/tests/test_flight.py
+++ b/tests/test_flight.py
@@ -604,12 +604,12 @@ def test_max_values(flight_calisto_robust):
     calculated by hand, it was just copied from the test results. This is
     because the expected values are not easy to calculate by hand, and the
     results are not expected to change. If the results change, the test will
-    fail, and the expected values must be updated. If if want to update the
-    values, always double check if the results are really correct. Acceptable
-    reasons for changes in the results are: 1) changes in the code that
-    improve the accuracy of the simulation, 2) a bug was found and fixed. Keep
-    in mind that other tests may be more accurate than this one, for example,
-    the acceptance tests, which are based on the results of real flights.
+    fail, and the expected values must be updated. If the values are updated,
+    always double check if the results are really correct. Acceptable reasons
+    for changes in the results are: 1) changes in the code that improve the
+    accuracy of the simulation, 2) a bug was found and fixed. Keep in mind that
+    other tests may be more accurate than this one, for example, the acceptance
+    tests, which are based on the results of real flights.
 
     Parameters
     ----------
@@ -622,7 +622,7 @@ def test_max_values(flight_calisto_robust):
     assert pytest.approx(105.2774, abs=atol) == test.max_acceleration_power_on
     assert pytest.approx(105.2774, abs=atol) == test.max_acceleration
     assert pytest.approx(0.85999, abs=atol) == test.max_mach_number
-    assert pytest.approx(285.90240, abs=atol) == test.max_speed
+    assert pytest.approx(285.94948, abs=atol) == test.max_speed
 
 
 def test_rail_buttons_forces(flight_calisto_custom_wind):
diff --git a/tests/unit/test_environment.py b/tests/unit/test_environment.py
index ac25533eb..8d676f426 100644
--- a/tests/unit/test_environment.py
+++ b/tests/unit/test_environment.py
@@ -1,10 +1,10 @@
+import json
 import os
 
 import numpy as np
 import numpy.ma as ma
 import pytest
 import pytz
-import json
 
 from rocketpy import Environment