diff --git a/.vscode/settings.json b/.vscode/settings.json
index 0af75e918..53e4c76d9 100644
--- a/.vscode/settings.json
+++ b/.vscode/settings.json
@@ -262,6 +262,7 @@
         "rtol",
         "rtype",
         "rucsoundings",
+        "runslow",
         "rwork",
         "savetxt",
         "savgol",
diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py
index 96afce730..71334e178 100644
--- a/rocketpy/environment/environment.py
+++ b/rocketpy/environment/environment.py
@@ -459,8 +459,7 @@ def __initialize_utm_coordinates(self):
                 flattening=self.ellipsoid.flattening,
                 semi_major_axis=self.ellipsoid.semi_major_axis,
             )
-        else:
-            # pragma: no cover
+        else:  # pragma: no cover
             warnings.warn(
                 "UTM coordinates are not available for latitudes "
                 "above 84 or below -80 degrees. The UTM conversions will fail."
@@ -715,8 +714,8 @@ def set_location(self, latitude, longitude):
 
         if not isinstance(latitude, NUMERICAL_TYPES) and isinstance(
             longitude, NUMERICAL_TYPES
-        ):
-            # pragma: no cover
+        ):  # pragma: no cover
+
             raise TypeError("Latitude and Longitude must be numbers!")
 
         # Store latitude and longitude
@@ -812,8 +811,8 @@ def max_expected_height(self):
 
     @max_expected_height.setter
     def max_expected_height(self, value):
-        if value < self.elevation:
-            raise ValueError(  # pragma: no cover
+        if value < self.elevation:  # pragma: no cover
+            raise ValueError(
                 "Max expected height cannot be lower than the surface elevation"
             )
         self._max_expected_height = value
@@ -952,8 +951,8 @@ def get_elevation_from_topographic_profile(self, lat, lon):
             Elevation provided by the topographic data, in meters.
         """
         # TODO: refactor this method.  pylint: disable=too-many-statements
-        if self.topographic_profile_activated is False:
-            raise ValueError(  # pragma: no cover
+        if self.topographic_profile_activated is False:  # pragma: no cover
+            raise ValueError(
                 "You must define a Topographic profile first, please use the "
                 "Environment.set_topographic_profile() method first."
             )
@@ -1285,8 +1284,8 @@ def set_atmospheric_model(  # pylint: disable=too-many-statements
                 self.process_forecast_reanalysis(dataset, dictionary)
             else:
                 self.process_ensemble(dataset, dictionary)
-        else:
-            raise ValueError(f"Unknown model type '{type}'.")  # pragma: no cover
+        else:  # pragma: no cover
+            raise ValueError(f"Unknown model type '{type}'.")
 
         if type not in ["ensemble"]:
             # Ensemble already computed these values
@@ -2578,7 +2577,7 @@ def set_earth_geometry(self, datum):
         }
         try:
             return ellipsoid[datum]
-        except KeyError as e:
+        except KeyError as e:  # pragma: no cover
             available_datums = ', '.join(ellipsoid.keys())
             raise AttributeError(
                 f"The reference system '{datum}' is not recognized. Please use one of "
@@ -2845,7 +2844,7 @@ def from_dict(cls, data):  # pylint: disable=too-many-statements
         return env
 
 
-if __name__ == "__main__":
+if __name__ == "__main__":  # pragma: no cover
     import doctest
 
     results = doctest.testmod()
diff --git a/rocketpy/environment/fetchers.py b/rocketpy/environment/fetchers.py
index 58537025b..b95f12118 100644
--- a/rocketpy/environment/fetchers.py
+++ b/rocketpy/environment/fetchers.py
@@ -79,12 +79,12 @@ def fetch_atmospheric_data_from_windy(lat, lon, model):
 
     try:
         response = requests.get(url).json()
-        if "data" not in response.keys():
+        if "data" not in response.keys():  # pragma: no cover
             raise ValueError(
                 f"Could not get a valid response for '{model}' from Windy. "
                 "Check if the coordinates are set inside the model's domain."
             )
-    except requests.exceptions.RequestException as e:
+    except requests.exceptions.RequestException as e:  # pragma: no cover
         if model == "iconEu":
             raise ValueError(
                 "Could not get a valid response for Icon-EU from Windy. "
@@ -315,8 +315,8 @@ def fetch_wyoming_sounding(file):
         If the response indicates the output format is invalid.
     """
     response = requests.get(file)
-    if response.status_code != 200:
-        raise ImportError(f"Unable to load {file}.")  # pragma: no cover
+    if response.status_code != 200:  # pragma: no cover
+        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]
@@ -330,7 +330,7 @@ def fetch_wyoming_sounding(file):
 
 
 @exponential_backoff(max_attempts=5, base_delay=2, max_delay=60)
-def fetch_noaaruc_sounding(file):
+def fetch_noaaruc_sounding(file):  # pragma: no cover
     """Fetches sounding data from a specified file using the NOAA RUC soundings.
 
     Parameters
diff --git a/rocketpy/environment/tools.py b/rocketpy/environment/tools.py
index 4fc3ca7c7..571a7de47 100644
--- a/rocketpy/environment/tools.py
+++ b/rocketpy/environment/tools.py
@@ -590,7 +590,7 @@ def utm_to_geodesic(  # pylint: disable=too-many-locals,too-many-statements
     return lat, lon
 
 
-if __name__ == "__main__":
+if __name__ == "__main__":  # pragma: no cover
     import doctest
 
     results = doctest.testmod()
diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py
index ccca0eab8..0b353c333 100644
--- a/rocketpy/mathutils/function.py
+++ b/rocketpy/mathutils/function.py
@@ -1481,7 +1481,7 @@ def plot(self, *args, **kwargs):
             else:
                 print("Error: Only functions with 1D or 2D domains can be plotted.")
 
-    def plot1D(self, *args, **kwargs):
+    def plot1D(self, *args, **kwargs):  # pragma: no cover
         """Deprecated method, use Function.plot_1d instead."""
         warnings.warn(
             "The `Function.plot1D` method is set to be deprecated and fully "
@@ -1581,7 +1581,7 @@ def plot_1d(  # pylint: disable=too-many-statements
         if return_object:
             return fig, ax
 
-    def plot2D(self, *args, **kwargs):
+    def plot2D(self, *args, **kwargs):  # pragma: no cover
         """Deprecated method, use Function.plot_2d instead."""
         warnings.warn(
             "The `Function.plot2D` method is set to be deprecated and fully "
@@ -2772,7 +2772,7 @@ def differentiate_complex_step(self, x, dx=1e-200, order=1):
         """
         if order == 1:
             return float(self.get_value_opt(x + dx * 1j).imag / dx)
-        else:
+        else:  # pragma: no cover
             raise NotImplementedError(
                 "Only 1st order derivatives are supported yet. Set order=1."
             )
@@ -3119,12 +3119,12 @@ def compose(self, func, extrapolate=False):
             The result of inputting the function into the function.
         """
         # Check if the input is a function
-        if not isinstance(func, Function):
+        if not isinstance(func, Function):  # pragma: no cover
             raise TypeError("Input must be a Function object.")
 
         if isinstance(self.source, np.ndarray) and isinstance(func.source, np.ndarray):
             # Perform bounds check for composition
-            if not extrapolate:
+            if not extrapolate:  # pragma: no cover
                 if func.min < self.x_initial or func.max > self.x_final:
                     raise ValueError(
                         f"Input Function image {func.min, func.max} must be within "
@@ -3197,7 +3197,7 @@ def savetxt(
 
         # create the datapoints
         if callable(self.source):
-            if lower is None or upper is None or samples is None:
+            if lower is None or upper is None or samples is None:  # pragma: no cover
                 raise ValueError(
                     "If the source is a callable, lower, upper and samples"
                     + " must be provided."
@@ -3264,7 +3264,7 @@ def __validate_source(self, source):  # pylint: disable=too-many-statements
                         self.__inputs__ = header[:-1]
                     if self.__outputs__ is None:
                         self.__outputs__ = [header[-1]]
-            except Exception as e:
+            except Exception as e:  # pragma: no cover
                 raise ValueError(
                     "Could not read the csv or txt file to create Function source."
                 ) from e
@@ -3323,6 +3323,7 @@ def __validate_inputs(self, inputs):
             if isinstance(inputs, (list, tuple)):
                 if len(inputs) == 1:
                     return inputs
+            # pragma: no cover
             raise ValueError(
                 "Inputs must be a string or a list of strings with "
                 "the length of the domain dimension."
@@ -3335,6 +3336,7 @@ def __validate_inputs(self, inputs):
                     isinstance(i, str) for i in inputs
                 ):
                     return inputs
+            # pragma: no cover
             raise ValueError(
                 "Inputs must be a list of strings with "
                 "the length of the domain dimension."
@@ -3611,7 +3613,7 @@ def reset_funcified_methods(instance):
             instance.__dict__.pop(key)
 
 
-if __name__ == "__main__":
+if __name__ == "__main__":  # pragma: no cover
     import doctest
 
     results = doctest.testmod()
diff --git a/rocketpy/mathutils/vector_matrix.py b/rocketpy/mathutils/vector_matrix.py
index 8147aa52e..ceac9a08b 100644
--- a/rocketpy/mathutils/vector_matrix.py
+++ b/rocketpy/mathutils/vector_matrix.py
@@ -1103,7 +1103,7 @@ def from_dict(cls, data):
         return cls(data)
 
 
-if __name__ == "__main__":
+if __name__ == "__main__":  # pragma: no cover
     import doctest
 
     results = doctest.testmod()
diff --git a/rocketpy/motors/fluid.py b/rocketpy/motors/fluid.py
index 62f91583b..e027702d0 100644
--- a/rocketpy/motors/fluid.py
+++ b/rocketpy/motors/fluid.py
@@ -30,9 +30,9 @@ def __post_init__(self):
             If the density is not a positive number.
         """
 
-        if not isinstance(self.name, str):
+        if not isinstance(self.name, str):  # pragma: no cover
             raise ValueError("The name must be a string.")
-        if self.density < 0:
+        if self.density < 0:  # pragma: no cover
             raise ValueError("The density must be a positive number.")
 
         # Initialize plots and prints object
diff --git a/rocketpy/motors/motor.py b/rocketpy/motors/motor.py
index ce72232f4..77c3a8633 100644
--- a/rocketpy/motors/motor.py
+++ b/rocketpy/motors/motor.py
@@ -247,7 +247,7 @@ class Function. Thrust units are Newtons.
             self._csys = 1
         elif coordinate_system_orientation == "combustion_chamber_to_nozzle":
             self._csys = -1
-        else:
+        else:  # pragma: no cover
             raise ValueError(
                 "Invalid coordinate system orientation. Options are "
                 "'nozzle_to_combustion_chamber' and 'combustion_chamber_to_nozzle'."
@@ -346,7 +346,7 @@ def burn_time(self, burn_time):
         else:
             if not callable(self.thrust.source):
                 self._burn_time = (self.thrust.x_array[0], self.thrust.x_array[-1])
-            else:
+            else:  # pragma: no cover
                 raise ValueError(
                     "When using a float or callable as thrust source, a burn_time"
                     " argument must be specified."
diff --git a/rocketpy/plots/rocket_plots.py b/rocketpy/plots/rocket_plots.py
index f305fcb99..6b571759b 100644
--- a/rocketpy/plots/rocket_plots.py
+++ b/rocketpy/plots/rocket_plots.py
@@ -385,7 +385,7 @@ def _draw_generic_surface(
             x_pos = position[2]
             # y position of the surface is the y position in the plot
             y_pos = position[1]
-        else:
+        else:  # pragma: no cover
             raise ValueError("Plane must be 'xz' or 'yz'.")
 
         ax.scatter(
@@ -633,7 +633,7 @@ def _draw_sensors(self, ax, sensors, plane):
                 # y position of the sensor is the y position in the plot
                 y_pos = pos[1]
                 normal_y = sensor.normal_vector.y
-            else:
+            else:  # pragma: no cover
                 raise ValueError("Plane must be 'xz' or 'yz'.")
 
             # line length is 2/5 of the rocket radius
diff --git a/rocketpy/prints/compare_prints.py b/rocketpy/prints/compare_prints.py
index 9b94f2da4..548dd0238 100644
--- a/rocketpy/prints/compare_prints.py
+++ b/rocketpy/prints/compare_prints.py
@@ -1,3 +1,3 @@
-class _ComparePrints:
+class _ComparePrints:  # pragma: no cover
     def __init__(self) -> None:
         pass
diff --git a/rocketpy/rocket/aero_surface/nose_cone.py b/rocketpy/rocket/aero_surface/nose_cone.py
index 7881b402c..954ce1ef8 100644
--- a/rocketpy/rocket/aero_surface/nose_cone.py
+++ b/rocketpy/rocket/aero_surface/nose_cone.py
@@ -135,7 +135,7 @@ def __init__(  # pylint: disable=too-many-statements
         self._base_radius = base_radius
         self._length = length
         if bluffness is not None:
-            if bluffness > 1 or bluffness < 0:
+            if bluffness > 1 or bluffness < 0:  # pragma: no cover
                 raise ValueError(
                     f"Bluffness ratio of {bluffness} is out of range. "
                     "It must be between 0 and 1."
@@ -286,7 +286,7 @@ def theta(x):
             self.y_nosecone = Function(
                 lambda x: self.base_radius * np.power(x / self.length, self.power)
             )
-        else:
+        else:  # pragma: no cover
             raise ValueError(
                 f"Nose Cone kind '{self.kind}' not found, "
                 + "please use one of the following Nose Cone kinds:"
@@ -317,12 +317,11 @@ def bluffness(self, value):
                 raise ValueError(
                     "Parameter 'bluffness' must be None or 0 when using a nose cone kind 'powerseries'."
                 )
-        if value is not None:
-            if value > 1 or value < 0:
-                raise ValueError(
-                    f"Bluffness ratio of {value} is out of range. "
-                    "It must be between 0 and 1."
-                )
+        if value is not None and not 0 <= value <= 1:  # pragma: no cover
+            raise ValueError(
+                f"Bluffness ratio of {value} is out of range. "
+                "It must be between 0 and 1."
+            )
         self._bluffness = value
         self.evaluate_nose_shape()
 
diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py
index c13714576..ef69f7641 100644
--- a/rocketpy/rocket/rocket.py
+++ b/rocketpy/rocket/rocket.py
@@ -278,7 +278,7 @@ def __init__(  # pylint: disable=too-many-statements
             self._csys = 1
         elif coordinate_system_orientation == "nose_to_tail":
             self._csys = -1
-        else:
+        else:  # pragma: no cover
             raise TypeError(
                 "Invalid coordinate system orientation. Please choose between "
                 + '"tail_to_nose" and "nose_to_tail".'
@@ -1173,7 +1173,7 @@ def add_nose(
         self.add_surfaces(nose, position)
         return nose
 
-    def add_fins(self, *args, **kwargs):
+    def add_fins(self, *args, **kwargs):  # pragma: no cover
         """See Rocket.add_trapezoidal_fins for documentation.
         This method is set to be deprecated in version 1.0.0 and fully removed
         by version 2.0.0. Use Rocket.add_trapezoidal_fins instead. It keeps the
diff --git a/rocketpy/sensitivity/sensitivity_model.py b/rocketpy/sensitivity/sensitivity_model.py
index de72cd0d2..428897bff 100644
--- a/rocketpy/sensitivity/sensitivity_model.py
+++ b/rocketpy/sensitivity/sensitivity_model.py
@@ -140,10 +140,6 @@ def set_target_variables_nominal(self, target_variables_nominal_value):
             self.target_variables_info[target_variable]["nominal_value"] = (
                 target_variables_nominal_value[i]
             )
-        for i, target_variable in enumerate(self.target_variables_names):
-            self.target_variables_info[target_variable]["nominal_value"] = (
-                target_variables_nominal_value[i]
-            )
 
         self._nominal_target_passed = True
 
@@ -356,12 +352,12 @@ def __check_requirements(self):
             version = ">=0" if not version else version
             try:
                 check_requirement_version(module_name, version)
-            except (ValueError, ImportError) as e:
+            except (ValueError, ImportError) as e:  # pragma: no cover
                 has_error = True
                 print(
                     f"The following error occurred while importing {module_name}: {e}"
                 )
-        if has_error:
+        if has_error:  # pragma: no cover
             print(
                 "Given the above errors, some methods may not work. Please run "
                 + "'pip install rocketpy[sensitivity]' to install extra requirements."
diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index a407e6928..941f4cf70 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -615,7 +615,7 @@ def __init__(  # pylint: disable=too-many-arguments,too-many-statements
         self.env = environment
         self.rocket = rocket
         self.rail_length = rail_length
-        if self.rail_length <= 0:
+        if self.rail_length <= 0:  # pragma: no cover
             raise ValueError("Rail length must be a positive value.")
         self.parachutes = self.rocket.parachutes[:]
         self.inclination = inclination
@@ -872,11 +872,11 @@ def __simulate(self, verbose):
                             for t_root in t_roots
                             if 0 < t_root.real < t1 and abs(t_root.imag) < 0.001
                         ]
-                        if len(valid_t_root) > 1:
+                        if len(valid_t_root) > 1:  # pragma: no cover
                             raise ValueError(
                                 "Multiple roots found when solving for rail exit time."
                             )
-                        if len(valid_t_root) == 0:
+                        if len(valid_t_root) == 0:  # pragma: no cover
                             raise ValueError(
                                 "No valid roots found when solving for rail exit time."
                             )
@@ -951,7 +951,7 @@ def __simulate(self, verbose):
                             for t_root in t_roots
                             if abs(t_root.imag) < 0.001 and 0 < t_root.real < t1
                         ]
-                        if len(valid_t_root) > 1:
+                        if len(valid_t_root) > 1:  # pragma: no cover
                             raise ValueError(
                                 "Multiple roots found when solving for impact time."
                             )
@@ -1226,7 +1226,7 @@ def __init_controllers(self):
         self._controllers = self.rocket._controllers[:]
         self.sensors = self.rocket.sensors.get_components()
         if self._controllers or self.sensors:
-            if self.time_overshoot:
+            if self.time_overshoot:  # pragma: no cover
                 self.time_overshoot = False
                 warnings.warn(
                     "time_overshoot has been set to False due to the presence "
@@ -1266,7 +1266,7 @@ def __set_ode_solver(self, solver):
         else:
             try:
                 self._solver = ODE_SOLVER_MAP[solver]
-            except KeyError as e:
+            except KeyError as e:  # pragma: no cover
                 raise ValueError(
                     f"Invalid ``ode_solver`` input: {solver}. "
                     f"Available options are: {', '.join(ODE_SOLVER_MAP.keys())}"
@@ -1398,7 +1398,7 @@ def udot_rail1(self, t, u, post_processing=False):
 
         return [vx, vy, vz, ax, ay, az, 0, 0, 0, 0, 0, 0, 0]
 
-    def udot_rail2(self, t, u, post_processing=False):
+    def udot_rail2(self, t, u, post_processing=False):  # pragma: no cover
         """[Still not implemented] Calculates derivative of u state vector with
         respect to time when rocket is flying in 3 DOF motion in the rail.
 
@@ -3531,7 +3531,7 @@ def __len__(self):
         def __repr__(self):
             return str(self.list)
 
-        def display_warning(self, *messages):
+        def display_warning(self, *messages):  # pragma: no cover
             """A simple function to print a warning message."""
             print("WARNING:", *messages)
 
diff --git a/rocketpy/tools.py b/rocketpy/tools.py
index a59cc1bd0..9962e9442 100644
--- a/rocketpy/tools.py
+++ b/rocketpy/tools.py
@@ -29,11 +29,10 @@
 
 
 def tuple_handler(value):
-    """Transforms the input value into a tuple that
-    represents a range. If the input is an int or float,
-    the output is a tuple from zero to the input value. If
-    the input is a tuple or list, the output is a tuple with
-    the same range.
+    """Transforms the input value into a tuple that represents a range. If the
+    input is an int or float, the output is a tuple from zero to the input
+    value. If the input is a tuple or list, the output is a tuple with the same
+    range.
 
     Parameters
     ----------
@@ -265,7 +264,7 @@ def get_distribution(distribution_function_name):
     }
     try:
         return distributions[distribution_function_name]
-    except KeyError as e:
+    except KeyError as e:  # pragma: no cover
         raise ValueError(
             f"Distribution function '{distribution_function_name}' not found, "
             + "please use one of the following np.random distribution function:"
@@ -915,7 +914,7 @@ def import_optional_dependency(name):
     """
     try:
         module = importlib.import_module(name)
-    except ImportError as exc:
+    except ImportError as exc:  # pragma: no cover
         module_name = name.split(".")[0]
         package_name = INSTALL_MAPPING.get(module_name, module_name)
         raise ImportError(
@@ -979,7 +978,8 @@ def wrapper(*args, **kwargs):
             for i in range(max_attempts):
                 try:
                     return func(*args, **kwargs)
-                except Exception as e:  # pylint: disable=broad-except
+                # pylint: disable=broad-except
+                except Exception as e:  # pragma: no cover
                     if i == max_attempts - 1:
                         raise e from None
                     delay = min(delay * 2, max_delay)
@@ -1205,7 +1205,7 @@ def from_hex_decode(obj_bytes, decoder=base64.b85decode):
     return dill.loads(decoder(bytes.fromhex(obj_bytes)))
 
 
-if __name__ == "__main__":
+if __name__ == "__main__":  # pragma: no cover
     import doctest
 
     res = doctest.testmod()
diff --git a/rocketpy/utilities.py b/rocketpy/utilities.py
index f6e68d162..4dcf7ebf5 100644
--- a/rocketpy/utilities.py
+++ b/rocketpy/utilities.py
@@ -14,7 +14,6 @@
 from .simulation.flight import Flight
 
 
-# TODO: Needs tests
 def compute_cd_s_from_drop_test(
     terminal_velocity, rocket_mass, air_density=1.225, g=9.80665
 ):
@@ -39,13 +38,34 @@ def compute_cd_s_from_drop_test(
     -------
     cd_s : float
         Number equal to drag coefficient times reference area for parachute.
+    """
+    return 2 * rocket_mass * g / ((terminal_velocity**2) * air_density)
+
 
+def check_constant(f, eps):
     """
+    Check for three consecutive elements in the list that are approximately
+    equal within a tolerance.
 
-    return 2 * rocket_mass * g / ((terminal_velocity**2) * air_density)
+    Parameters
+    ----------
+    f : list or array
+        A list or array of numerical values.
+    eps : float
+        The tolerance level for comparing the elements.
+
+    Returns
+    -------
+    int or None
+        The index of the first element in the first sequence of three
+        consecutive elements that are approximately equal within the tolerance.
+        Returns None if no such sequence is found.
+    """
+    for i in range(len(f) - 2):
+        if abs(f[i + 2] - f[i + 1]) < eps and abs(f[i + 1] - f[i]) < eps:
+            return i
 
 
-# TODO: Needs tests
 def calculate_equilibrium_altitude(
     rocket_mass,
     cd_s,
@@ -90,7 +110,6 @@ def calculate_equilibrium_altitude(
         affect the final result if the value is not high enough. Increase the
         estimative in case the final solution is not founded.
 
-
     Returns
     -------
     altitude_function: Function
@@ -103,30 +122,8 @@ def calculate_equilibrium_altitude(
     """
     final_sol = {}
 
-    if v0 >= 0:
-        print("Please set a valid negative value for v0")
-        return None
-
-    # TODO: Improve docs
-    def check_constant(f, eps):
-        """_summary_
-
-        Parameters
-        ----------
-        f : array, list
-
-            _description_
-        eps : float
-            _description_
-
-        Returns
-        -------
-        int, None
-            _description_
-        """
-        for i in range(len(f) - 2):
-            if abs(f[i + 2] - f[i + 1]) < eps and abs(f[i + 1] - f[i]) < eps:
-                return i
+    if v0 >= 0:  # pragma: no cover
+        raise ValueError("Please set a valid negative value for v0")
 
     if env is None:
         environment = Environment(
@@ -138,21 +135,20 @@ def check_constant(f, eps):
     else:
         environment = env
 
-    # TODO: Improve docs
     def du(z, u):
-        """_summary_
+        """Returns the derivative of the velocity at a given altitude.
 
         Parameters
         ----------
         z : float
-            _description_
+            altitude, in meters, at a given time
         u : float
             velocity, in m/s, at a given z altitude
 
         Returns
         -------
         float
-            _description_
+            velocity at a given altitude
         """
         return (
             u[1],
@@ -258,7 +254,7 @@ def fin_flutter_analysis(
                 found_fin = True
             else:
                 warnings.warn("More than one fin set found. The last one will be used.")
-    if not found_fin:
+    if not found_fin:  # pragma: no cover
         raise AttributeError(
             "There is no TrapezoidalFins in the rocket, can't run Flutter Analysis."
         )
@@ -442,7 +438,7 @@ def _flutter_prints(
 
 
 # TODO: deprecate and delete this function. Never used and now we have Monte Carlo.
-def create_dispersion_dictionary(filename):
+def create_dispersion_dictionary(filename):  # pragma: no cover
     """Creates a dictionary with the rocket data provided by a .csv file.
     File should be organized in four columns: attribute_class, parameter_name,
     mean_value, standard_deviation. The first row should be the header.
diff --git a/tests/fixtures/flight/flight_fixtures.py b/tests/fixtures/flight/flight_fixtures.py
index 0a9a39ff9..44f8b7c09 100644
--- a/tests/fixtures/flight/flight_fixtures.py
+++ b/tests/fixtures/flight/flight_fixtures.py
@@ -93,6 +93,22 @@ def flight_calisto_robust(calisto_robust, example_spaceport_env):
     )
 
 
+@pytest.fixture
+def flight_calisto_robust_solid_eom(calisto_robust, example_spaceport_env):
+    """Similar to flight_calisto_robust, but with the equations of motion set to
+    "solid_propulsion".
+    """
+    return Flight(
+        environment=example_spaceport_env,
+        rocket=calisto_robust,
+        rail_length=5.2,
+        inclination=85,
+        heading=0,
+        terminate_on_apogee=False,
+        equations_of_motion="solid_propulsion",
+    )
+
+
 @pytest.fixture
 def flight_calisto_liquid_modded(calisto_liquid_modded, example_plain_env):
     """A rocketpy.Flight object of the Calisto rocket modded for a liquid
diff --git a/tests/fixtures/surfaces/surface_fixtures.py b/tests/fixtures/surfaces/surface_fixtures.py
index 396206bd7..bf6e384c4 100644
--- a/tests/fixtures/surfaces/surface_fixtures.py
+++ b/tests/fixtures/surfaces/surface_fixtures.py
@@ -1,7 +1,13 @@
 import pytest
 
-from rocketpy import NoseCone, RailButtons, Tail, TrapezoidalFins
-from rocketpy.rocket.aero_surface.fins.free_form_fins import FreeFormFins
+from rocketpy.rocket.aero_surface import (
+    EllipticalFins,
+    FreeFormFins,
+    NoseCone,
+    RailButtons,
+    Tail,
+    TrapezoidalFins,
+)
 
 
 @pytest.fixture
@@ -94,3 +100,16 @@ def calisto_rail_buttons():
         angular_position=45,
         name="Rail Buttons",
     )
+
+
+@pytest.fixture
+def elliptical_fin_set():
+    return EllipticalFins(
+        n=4,
+        span=0.100,
+        root_chord=0.120,
+        rocket_radius=0.0635,
+        cant_angle=0,
+        airfoil=None,
+        name="Test Elliptical Fins",
+    )
diff --git a/tests/integration/test_flight.py b/tests/integration/test_flight.py
index 7edf60358..3efc0c285 100644
--- a/tests/integration/test_flight.py
+++ b/tests/integration/test_flight.py
@@ -10,9 +10,12 @@
 plt.rcParams.update({"figure.max_open_warning": 0})
 
 
+@pytest.mark.parametrize(
+    "flight_fixture", ["flight_calisto_robust", "flight_calisto_robust_solid_eom"]
+)
 @patch("matplotlib.pyplot.show")
 # pylint: disable=unused-argument
-def test_all_info(mock_show, flight_calisto_robust):
+def test_all_info(mock_show, request, flight_fixture):
     """Test that the flight class is working as intended. This basically calls
     the all_info() method and checks if it returns None. It is not testing if
     the values are correct, but whether the method is working without errors.
@@ -21,11 +24,13 @@ def test_all_info(mock_show, flight_calisto_robust):
     ----------
     mock_show : unittest.mock.MagicMock
         Mock object to replace matplotlib.pyplot.show
-    flight_calisto_robust : rocketpy.Flight
-        Flight object to be tested. See the conftest.py file for more info
-        regarding this pytest fixture.
+    request : _pytest.fixtures.FixtureRequest
+        Request object to access the fixture dynamically.
+    flight_fixture : str
+        Name of the flight fixture to be tested.
     """
-    assert flight_calisto_robust.all_info() is None
+    flight = request.getfixturevalue(flight_fixture)
+    assert flight.all_info() is None
 
 
 @pytest.mark.slow
diff --git a/tests/unit/test_aero_surfaces.py b/tests/unit/test_aero_surfaces.py
index 5258814db..a59820b5f 100644
--- a/tests/unit/test_aero_surfaces.py
+++ b/tests/unit/test_aero_surfaces.py
@@ -1,3 +1,5 @@
+from unittest.mock import patch
+
 import pytest
 
 from rocketpy import NoseCone
@@ -71,3 +73,69 @@ def test_powerseries_nosecones_setters(power, invalid_power, new_power):
     expected_k = (2 * new_power) / ((2 * new_power) + 1)
 
     assert pytest.approx(test_nosecone.k) == expected_k
+
+
+@patch("matplotlib.pyplot.show")
+def test_elliptical_fins_draw(
+    mock_show, elliptical_fin_set
+):  # pylint: disable=unused-argument
+    assert elliptical_fin_set.plots.draw(filename=None) is None
+
+
+def test_nose_cone_info(calisto_nose_cone):
+    assert calisto_nose_cone.info() is None
+
+
+@patch("matplotlib.pyplot.show")
+def test_nose_cone_draw(
+    mock_show, calisto_nose_cone
+):  # pylint: disable=unused-argument
+    assert calisto_nose_cone.draw(filename=None) is None
+
+
+def test_trapezoidal_fins_info(calisto_trapezoidal_fins):
+    assert calisto_trapezoidal_fins.info() is None
+
+
+def test_trapezoidal_fins_tip_chord_setter(calisto_trapezoidal_fins):
+    calisto_trapezoidal_fins.tip_chord = 0.1
+    assert calisto_trapezoidal_fins.tip_chord == 0.1
+
+
+def test_trapezoidal_fins_root_chord_setter(calisto_trapezoidal_fins):
+    calisto_trapezoidal_fins.root_chord = 0.1
+    assert calisto_trapezoidal_fins.root_chord == 0.1
+
+
+def test_trapezoidal_fins_sweep_angle_setter(calisto_trapezoidal_fins):
+    calisto_trapezoidal_fins.sweep_angle = 0.1
+    assert calisto_trapezoidal_fins.sweep_angle == 0.1
+
+
+def test_trapezoidal_fins_sweep_length_setter(calisto_trapezoidal_fins):
+    calisto_trapezoidal_fins.sweep_length = 0.1
+    assert calisto_trapezoidal_fins.sweep_length == 0.1
+
+
+def test_tail_info(calisto_tail):
+    assert calisto_tail.info() is None
+
+
+def test_tail_length_setter(calisto_tail):
+    calisto_tail.length = 0.1
+    assert calisto_tail.length == 0.1
+
+
+def test_tail_rocket_radius_setter(calisto_tail):
+    calisto_tail.rocket_radius = 0.1
+    assert calisto_tail.rocket_radius == 0.1
+
+
+def test_tail_bottom_radius_setter(calisto_tail):
+    calisto_tail.bottom_radius = 0.1
+    assert calisto_tail.bottom_radius == 0.1
+
+
+def test_tail_top_radius_setter(calisto_tail):
+    calisto_tail.top_radius = 0.1
+    assert calisto_tail.top_radius == 0.1
diff --git a/tests/unit/test_flight_time_nodes.py b/tests/unit/test_flight_time_nodes.py
index dcdc11eff..20769b1f8 100644
--- a/tests/unit/test_flight_time_nodes.py
+++ b/tests/unit/test_flight_time_nodes.py
@@ -99,3 +99,13 @@ def test_time_node_lt(flight_calisto):
     node2 = flight_calisto.TimeNodes.TimeNode(2.0, [], [], [])
     assert node1 < node2
     assert not node2 < node1
+
+
+def test_time_node_repr(flight_calisto):
+    node = flight_calisto.TimeNodes.TimeNode(1.0, [], [], [])
+    assert isinstance(repr(node), str)
+
+
+def test_time_nodes_repr(flight_calisto):
+    time_nodes = flight_calisto.TimeNodes()
+    assert isinstance(repr(time_nodes), str)
diff --git a/tests/unit/test_function.py b/tests/unit/test_function.py
index df540c1ae..e138862c0 100644
--- a/tests/unit/test_function.py
+++ b/tests/unit/test_function.py
@@ -787,3 +787,107 @@ def test_low_pass_filter(alpha):
             f"The filtered value at index {i} is not the expected value. "
             f"Expected: {expected}, Actual: {filtered_func.source[i][1]}"
         )
+
+
+def test_average_function_ndarray():
+
+    dummy_function = Function(
+        source=[
+            [0, 0],
+            [1, 1],
+            [2, 0],
+            [3, 1],
+            [4, 0],
+            [5, 1],
+            [6, 0],
+            [7, 1],
+            [8, 0],
+            [9, 1],
+        ],
+        inputs=["x"],
+        outputs=["y"],
+    )
+    avg_function = dummy_function.average_function()
+
+    assert isinstance(avg_function, Function)
+    assert np.isclose(avg_function(0), 0)
+    assert np.isclose(avg_function(9), 0.5)
+
+
+def test_average_function_callable():
+
+    dummy_function = Function(lambda x: 2)
+    avg_function = dummy_function.average_function(lower=0)
+
+    assert isinstance(avg_function, Function)
+    assert np.isclose(avg_function(1), 2)
+    assert np.isclose(avg_function(9), 2)
+
+
+@pytest.mark.parametrize(
+    "lower, upper, sampling_frequency, window_size, step_size, remove_dc, only_positive",
+    [
+        (0, 10, 100, 1, 0.5, True, True),
+        (0, 10, 100, 1, 0.5, True, False),
+        (0, 10, 100, 1, 0.5, False, True),
+        (0, 10, 100, 1, 0.5, False, False),
+        (0, 20, 200, 2, 1, True, True),
+    ],
+)
+def test_short_time_fft(
+    lower, upper, sampling_frequency, window_size, step_size, remove_dc, only_positive
+):
+    """Test the short_time_fft method of the Function class.
+
+    Parameters
+    ----------
+    lower : float
+        Lower bound of the time range.
+    upper : float
+        Upper bound of the time range.
+    sampling_frequency : float
+        Sampling frequency at which to perform the Fourier transform.
+    window_size : float
+        Size of the window for the STFT, in seconds.
+    step_size : float
+        Step size for the window, in seconds.
+    remove_dc : bool
+        If True, the DC component is removed from each window before
+        computing the Fourier transform.
+    only_positive: bool
+        If True, only the positive frequencies are returned.
+    """
+    # Generate a test signal
+    t = np.linspace(lower, upper, int((upper - lower) * sampling_frequency))
+    signal = np.sin(2 * np.pi * 5 * t)  # 5 Hz sine wave
+    func = Function(np.column_stack((t, signal)))
+
+    # Perform STFT
+    stft_results = func.short_time_fft(
+        lower=lower,
+        upper=upper,
+        sampling_frequency=sampling_frequency,
+        window_size=window_size,
+        step_size=step_size,
+        remove_dc=remove_dc,
+        only_positive=only_positive,
+    )
+
+    # Check the results
+    assert isinstance(stft_results, list)
+    assert all(isinstance(f, Function) for f in stft_results)
+
+    for f in stft_results:
+        assert f.get_inputs() == ["Frequency (Hz)"]
+        assert f.get_outputs() == ["Amplitude"]
+        assert f.get_interpolation_method() == "linear"
+        assert f.get_extrapolation_method() == "zero"
+
+        frequencies = f.source[:, 0]
+        # amplitudes = f.source[:, 1]
+
+        if only_positive:
+            assert np.all(frequencies >= 0)
+        else:
+            assert np.all(frequencies >= -sampling_frequency / 2)
+            assert np.all(frequencies <= sampling_frequency / 2)
diff --git a/tests/unit/test_sensitivity.py b/tests/unit/test_sensitivity.py
index 86f8a918c..af8c69ac7 100644
--- a/tests/unit/test_sensitivity.py
+++ b/tests/unit/test_sensitivity.py
@@ -1,13 +1,11 @@
+from unittest.mock import patch
+
 import numpy as np
 import pytest
 
 from rocketpy.sensitivity import SensitivityModel
 
-# TODO: for some weird reason, these tests are not passing in the CI, but
-#       passing locally. Need to investigate why.
-
 
-@pytest.mark.skip(reason="legacy test")
 def test_initialization():
     parameters_names = ["param1", "param2"]
     target_variables_names = ["target1", "target2"]
@@ -21,7 +19,6 @@ def test_initialization():
     assert not model._fitted
 
 
-@pytest.mark.skip(reason="legacy test")
 def test_set_parameters_nominal():
     parameters_names = ["param1", "param2"]
     target_variables_names = ["target1", "target2"]
@@ -35,8 +32,16 @@ def test_set_parameters_nominal():
     assert model.parameters_info["param1"]["nominal_mean"] == 1.0
     assert model.parameters_info["param2"]["nominal_sd"] == 0.2
 
+    # check dimensions mismatch error raise
+    incorrect_nominal_mean = np.array([1.0])
+    with pytest.raises(ValueError):
+        model.set_parameters_nominal(incorrect_nominal_mean, parameters_nominal_sd)
+
+    incorrect_nominal_sd = np.array([0.1])
+    with pytest.raises(ValueError):
+        model.set_parameters_nominal(parameters_nominal_mean, incorrect_nominal_sd)
+
 
-@pytest.mark.skip(reason="legacy test")
 def test_set_target_variables_nominal():
     parameters_names = ["param1", "param2"]
     target_variables_names = ["target1", "target2"]
@@ -49,9 +54,13 @@ def test_set_target_variables_nominal():
     assert model.target_variables_info["target1"]["nominal_value"] == 10.0
     assert model.target_variables_info["target2"]["nominal_value"] == 20.0
 
+    # check dimensions mismatch error raise
+    incorrect_nominal_value = np.array([10.0])
+    with pytest.raises(ValueError):
+        model.set_target_variables_nominal(incorrect_nominal_value)
+
 
-@pytest.mark.skip(reason="legacy test")
-def test_fit_method():
+def test_fit_method_one_target():
     parameters_names = ["param1", "param2"]
     target_variables_names = ["target1"]
     model = SensitivityModel(parameters_names, target_variables_names)
@@ -65,7 +74,20 @@ def test_fit_method():
     assert model.number_of_samples == 3
 
 
-@pytest.mark.skip(reason="legacy test")
+def test_fit_method_multiple_target():
+    parameters_names = ["param1", "param2"]
+    target_variables_names = ["target1", "target2"]
+    model = SensitivityModel(parameters_names, target_variables_names)
+
+    parameters_matrix = np.array([[1.0, 2.0], [2.0, 3.0], [3.0, 4.0]])
+    target_data = np.array([[10.0, 12.0, 14.0], [11.0, 13.0, 17.0]]).T
+
+    model.fit(parameters_matrix, target_data)
+
+    assert model._fitted
+    assert model.number_of_samples == 3
+
+
 def test_fit_raises_error_on_mismatched_dimensions():
     parameters_names = ["param1", "param2"]
     target_variables_names = ["target1"]
@@ -78,7 +100,6 @@ def test_fit_raises_error_on_mismatched_dimensions():
         model.fit(parameters_matrix, target_data)
 
 
-@pytest.mark.skip(reason="legacy test")
 def test_check_conformity():
     parameters_names = ["param1", "param2"]
     target_variables_names = ["target1", "target2"]
@@ -90,7 +111,6 @@ def test_check_conformity():
     model._SensitivityModel__check_conformity(parameters_matrix, target_data)
 
 
-@pytest.mark.skip(reason="legacy test")
 def test_check_conformity_raises_error():
     parameters_names = ["param1", "param2"]
     target_variables_names = ["target1", "target2"]
@@ -101,3 +121,42 @@ def test_check_conformity_raises_error():
 
     with pytest.raises(ValueError):
         model._SensitivityModel__check_conformity(parameters_matrix, target_data)
+
+    parameters_matrix2 = np.array([[1.0, 2.0, 3.0], [2.0, 3.0, 4.0]])
+
+    with pytest.raises(ValueError):
+        model._SensitivityModel__check_conformity(parameters_matrix2, target_data)
+
+    target_data2 = np.array([10.0, 12.0])
+
+    with pytest.raises(ValueError):
+        model._SensitivityModel__check_conformity(parameters_matrix, target_data2)
+
+    target_variables_names = ["target1"]
+    model = SensitivityModel(parameters_names, target_variables_names)
+
+    target_data = np.array([[10.0, 20.0], [12.0, 22.0], [14.0, 24.0]])
+
+    with pytest.raises(ValueError):
+        model._SensitivityModel__check_conformity(parameters_matrix, target_data)
+
+
+@patch("matplotlib.pyplot.show")
+def test_prints_and_plots(mock_show):  # pylint: disable=unused-argument
+    parameters_names = ["param1", "param2"]
+    target_variables_names = ["target1"]
+    model = SensitivityModel(parameters_names, target_variables_names)
+
+    parameters_matrix = np.array([[1.0, 2.0], [2.0, 3.0], [3.0, 4.0]])
+    target_data = np.array([10.0, 12.0, 14.0])
+
+    # tests if an error is raised if summary is called before print
+    with pytest.raises(ValueError):
+        model.info()
+
+    model.fit(parameters_matrix, target_data)
+    assert model.all_info() is None
+
+    nominal_target = np.array([12.0])
+    model.set_target_variables_nominal(nominal_target)
+    assert model.all_info() is None
diff --git a/tests/unit/test_tank.py b/tests/unit/test_tank.py
index a0152cdfe..95179ccfa 100644
--- a/tests/unit/test_tank.py
+++ b/tests/unit/test_tank.py
@@ -1,10 +1,13 @@
 from math import isclose
 from pathlib import Path
+from unittest.mock import patch
 
 import numpy as np
 import pytest
 import scipy.integrate as spi
 
+from rocketpy.motors import TankGeometry
+
 BASE_PATH = Path("./data/rockets/berkeley/")
 
 
@@ -355,3 +358,8 @@ def expected_gas_inertia(t):
         atol=1e-3,
         rtol=1e-2,
     )
+
+
+@patch("matplotlib.pyplot.show")
+def test_tank_geometry_plots_info(mock_show):  # pylint: disable=unused-argument
+    assert TankGeometry({(0, 5): 1}).plots.all() is None
diff --git a/tests/unit/test_tools.py b/tests/unit/test_tools.py
index 9b321ea0e..fcf67ad37 100644
--- a/tests/unit/test_tools.py
+++ b/tests/unit/test_tools.py
@@ -6,6 +6,7 @@
     euler313_to_quaternions,
     find_roots_cubic_function,
     haversine,
+    tuple_handler,
 )
 
 
@@ -72,3 +73,30 @@ def test_cardanos_root_finding():
 def test_haversine(lat0, lon0, lat1, lon1, expected_distance):
     distance = haversine(lat0, lon0, lat1, lon1)
     assert np.isclose(distance, expected_distance, rtol=1e-2)
+
+
+@pytest.mark.parametrize(
+    "input_value, expected_output",
+    [
+        (5, (0, 5)),
+        (3.5, (0, 3.5)),
+        ([7], (0, 7)),
+        ((8,), (0, 8)),
+        ([2, 4], (2, 4)),
+        ((1, 3), (1, 3)),
+    ],
+)
+def test_tuple_handler(input_value, expected_output):
+    assert tuple_handler(input_value) == expected_output
+
+
+@pytest.mark.parametrize(
+    "input_value, expected_exception",
+    [
+        ([1, 2, 3], ValueError),
+        ((4, 5, 6), ValueError),
+    ],
+)
+def test_tuple_handler_exceptions(input_value, expected_exception):
+    with pytest.raises(expected_exception):
+        tuple_handler(input_value)
diff --git a/tests/unit/test_units.py b/tests/unit/test_units.py
new file mode 100644
index 000000000..76cb7ed89
--- /dev/null
+++ b/tests/unit/test_units.py
@@ -0,0 +1,84 @@
+import pytest
+
+from rocketpy.units import conversion_factor, convert_temperature, convert_units
+
+
+class TestConvertTemperature:
+    """Tests for the convert_temperature function."""
+
+    def test_convert_temperature_same_unit(self):
+        assert convert_temperature(300, "K", "K") == 300
+        assert convert_temperature(27, "degC", "degC") == 27
+        assert convert_temperature(80, "degF", "degF") == 80
+
+    def test_convert_temperature_kelvin_to_celsius(self):
+        assert convert_temperature(300, "K", "degC") == pytest.approx(26.85, rel=1e-2)
+
+    def test_convert_temperature_kelvin_to_fahrenheit(self):
+        assert convert_temperature(300, "K", "degF") == pytest.approx(80.33, rel=1e-2)
+
+    def test_convert_temperature_celsius_to_kelvin(self):
+        assert convert_temperature(27, "degC", "K") == pytest.approx(300.15, rel=1e-2)
+
+    def test_convert_temperature_celsius_to_fahrenheit(self):
+        assert convert_temperature(27, "degC", "degF") == pytest.approx(80.6, rel=1e-2)
+
+    def test_convert_temperature_fahrenheit_to_kelvin(self):
+        assert convert_temperature(80, "degF", "K") == pytest.approx(299.817, rel=1e-2)
+
+    def test_convert_temperature_fahrenheit_to_celsius(self):
+        assert convert_temperature(80, "degF", "degC") == pytest.approx(26.67, rel=1e-2)
+
+    def test_convert_temperature_invalid_conversion(self):
+        with pytest.raises(ValueError):
+            convert_temperature(300, "K", "invalid_unit")
+        with pytest.raises(ValueError):
+            convert_temperature(300, "invalid_unit", "K")
+
+
+class TestConversionFactor:
+    """Tests for the conversion_factor function."""
+
+    def test_conversion_factor_same_unit(self):
+        assert conversion_factor("m", "m") == 1
+        assert conversion_factor("ft", "ft") == 1
+        assert conversion_factor("s", "s") == 1
+
+    def test_conversion_factor_m_to_ft(self):
+        assert conversion_factor("m", "ft") == pytest.approx(3.28084, rel=1e-2)
+
+    def test_conversion_factor_ft_to_m(self):
+        assert conversion_factor("ft", "m") == pytest.approx(0.3048, rel=1e-2)
+
+    def test_conversion_factor_s_to_min(self):
+        assert conversion_factor("s", "min") == pytest.approx(1 / 60, rel=1e-2)
+
+    def test_conversion_factor_min_to_s(self):
+        assert conversion_factor("min", "s") == pytest.approx(60, rel=1e-2)
+
+    def test_conversion_factor_invalid_conversion(self):
+        with pytest.raises(ValueError):
+            conversion_factor("m", "invalid_unit")
+        with pytest.raises(ValueError):
+            conversion_factor("invalid_unit", "m")
+
+
+class TestConvertUnits:
+    """Tests for the convert_units function."""
+
+    def test_convert_units_same_unit(self):
+        assert convert_units(300, "K", "K") == 300
+        assert convert_units(27, "degC", "degC") == 27
+        assert convert_units(80, "degF", "degF") == 80
+
+    def test_convert_units_kelvin_to_celsius(self):
+        assert convert_units(300, "K", "degC") == pytest.approx(26.85, rel=1e-2)
+
+    def test_convert_units_kelvin_to_fahrenheit(self):
+        assert convert_units(300, "K", "degF") == pytest.approx(80.33, rel=1e-2)
+
+    def test_convert_units_kilogram_to_pound(self):
+        assert convert_units(1, "kg", "lb") == pytest.approx(2.20462, rel=1e-2)
+
+    def test_convert_units_kilometer_to_mile(self):
+        assert convert_units(1, "km", "mi") == pytest.approx(0.621371, rel=1e-2)
diff --git a/tests/unit/test_utilities.py b/tests/unit/test_utilities.py
index a6d1972a7..67cff10e3 100644
--- a/tests/unit/test_utilities.py
+++ b/tests/unit/test_utilities.py
@@ -178,3 +178,34 @@ def test_get_instance_attributes(flight_calisto_robust):
             assert np.allclose(attr, value)
         else:
             assert attr == value
+
+
+@pytest.mark.parametrize(
+    "f, eps, expected",
+    [
+        ([1.0, 1.0, 1.0, 2.0, 3.0], 1e-6, 0),
+        ([1.0, 1.0, 1.0, 2.0, 3.0], 1e-1, 0),
+        ([1.0, 1.1, 1.2, 2.0, 3.0], 1e-1, None),
+        ([1.0, 1.0, 1.0, 1.0, 1.0], 1e-6, 0),
+        ([1.0, 1.0, 1.0, 1.0, 1.0], 1e-1, 0),
+        ([1.0, 1.0, 1.0], 1e-6, 0),
+        ([1.0, 1.0], 1e-6, None),
+        ([1.0], 1e-6, None),
+        ([], 1e-6, None),
+    ],
+)
+def test_check_constant(f, eps, expected):
+    """Test if the function `check_constant` returns the correct index or None
+    for different scenarios.
+
+    Parameters
+    ----------
+    f : list or array
+        A list or array of numerical values.
+    eps : float
+        The tolerance level for comparing the elements.
+    expected : int or None
+        The expected result of the function.
+    """
+    result = utilities.check_constant(f, eps)
+    assert result == expected