diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1137336a0..c8e22a07a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -32,6 +32,7 @@ straightforward as possible.
 
 ### Added
 
+- ENH: Function Support for CSV Header Inputs [#542](https://github.com/RocketPy-Team/RocketPy/pull/542)
 - ENH: Shepard Optimized Interpolation - Multiple Inputs Support [#515](https://github.com/RocketPy-Team/RocketPy/pull/515)
 - ENH: adds new Function.savetxt method [#514](https://github.com/RocketPy-Team/RocketPy/pull/514)
 - ENH: Argument for Optional Mutation on Function Discretize [#519](https://github.com/RocketPy-Team/RocketPy/pull/519)
diff --git a/rocketpy/environment/environment_analysis.py b/rocketpy/environment/environment_analysis.py
index 2be707d73..c15b32551 100644
--- a/rocketpy/environment/environment_analysis.py
+++ b/rocketpy/environment/environment_analysis.py
@@ -901,10 +901,10 @@ def __parse_surface_data(self):
             # Extract data from weather file
             indices = (time_index, lon_index, lat_index)
             for key, value in surface_file_dict.items():
-                dictionary[date_string][hour_string][
-                    key
-                ] = self.__extract_surface_data_value(
-                    surface_data, value, indices, lon_array, lat_array
+                dictionary[date_string][hour_string][key] = (
+                    self.__extract_surface_data_value(
+                        surface_data, value, indices, lon_array, lat_array
+                    )
                 )
 
         # Get elevation, time index does not matter, use last one
diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py
index 7981bbbfa..048b125e2 100644
--- a/rocketpy/mathutils/function.py
+++ b/rocketpy/mathutils/function.py
@@ -3,6 +3,7 @@
 and more. This is a core class of our package, and should be maintained
 carefully as it may impact all the rest of the project.
 """
+
 import warnings
 from collections.abc import Iterable
 from copy import deepcopy
@@ -58,7 +59,8 @@ def __init__(
             and 'z' is the output.
 
             - string: Path to a CSV file. The file is read and converted into an
-            ndarray. The file can optionally contain a single header line.
+            ndarray. The file can optionally contain a single header line, see
+            notes below for more information.
 
             - Function: Copies the source of the provided Function object,
             creating a new Function with adjusted inputs and outputs.
@@ -93,12 +95,19 @@ def __init__(
 
         Notes
         -----
-        (I) CSV files can optionally contain a single header line. If present,
-        the header is ignored during processing.
-        (II) Fields in CSV files may be enclosed in double quotes. If fields are
-        not quoted, double quotes should not appear inside them.
+        (I) CSV files may include an optional single header line. If this
+        header line is present and contains names for each data column, those
+        names will be used to label the inputs and outputs unless specified
+        otherwise by the `inputs` and `outputs` arguments.
+        If the header is specified for only a few columns, it is ignored.
+
+        Commas in a header will be interpreted as a delimiter, which may cause
+        undesired input or output labeling. To avoid this, specify each input
+        and output name using the `inputs` and `outputs` arguments.
+
+        (II) Fields in CSV files may be enclosed in double quotes. If fields
+        are not quoted, double quotes should not appear inside them.
         """
-        # Set input and output
         if inputs is None:
             inputs = ["Scalar"]
         if outputs is None:
@@ -183,10 +192,18 @@ def set_source(self, source):
 
         Notes
         -----
-        (I) CSV files can optionally contain a single header line. If present,
-        the header is ignored during processing.
-        (II) Fields in CSV files may be enclosed in double quotes. If fields are
-        not quoted, double quotes should not appear inside them.
+        (I) CSV files may include an optional single header line. If this
+        header line is present and contains names for each data column, those
+        names will be used to label the inputs and outputs unless specified
+        otherwise. If the header is specified for only a few columns, it is
+        ignored.
+
+        Commas in a header will be interpreted as a delimiter, which may cause
+        undesired input or output labeling. To avoid this, specify each input
+        and output name using the `inputs` and `outputs` arguments.
+
+        (II) Fields in CSV files may be enclosed in double quotes. If fields
+        are not quoted, double quotes should not appear inside them.
 
         Returns
         -------
@@ -3018,7 +3035,7 @@ def _check_user_input(
         if isinstance(inputs, str):
             inputs = [inputs]
 
-        elif len(outputs) > 1:
+        if len(outputs) > 1:
             raise ValueError(
                 "Output must either be a string or have dimension 1, "
                 + f"it currently has dimension ({len(outputs)})."
@@ -3035,8 +3052,19 @@ def _check_user_input(
                 try:
                     source = np.loadtxt(source, delimiter=",", dtype=float)
                 except ValueError:
-                    # Skip header
-                    source = np.loadtxt(source, delimiter=",", dtype=float, skiprows=1)
+                    with open(source, "r") as file:
+                        header, *data = file.read().splitlines()
+
+                    header = [
+                        label.strip("'").strip('"') for label in header.split(",")
+                    ]
+                    source = np.loadtxt(data, delimiter=",", dtype=float)
+
+                    if len(source[0]) == len(header):
+                        if inputs == ["Scalar"]:
+                            inputs = header[:-1]
+                        if outputs == ["Scalar"]:
+                            outputs = [header[-1]]
                 except Exception as e:
                     raise ValueError(
                         "The source file is not a valid csv or txt file."
@@ -3054,7 +3082,7 @@ def _check_user_input(
 
             ## single dimension
             if source_dim == 2:
-                # possible interpolation values: llinear, polynomial, akima and spline
+                # possible interpolation values: linear, polynomial, akima and spline
                 if interpolation is None:
                     interpolation = "spline"
                 elif interpolation.lower() not in [
@@ -3105,7 +3133,7 @@ def _check_user_input(
             in_out_dim = len(inputs) + len(outputs)
             if source_dim != in_out_dim:
                 raise ValueError(
-                    "Source dimension ({source_dim}) does not match input "
+                    f"Source dimension ({source_dim}) does not match input "
                     + f"and output dimension ({in_out_dim})."
                 )
         return inputs, outputs, interpolation, extrapolation
diff --git a/rocketpy/plots/flight_plots.py b/rocketpy/plots/flight_plots.py
index e3b144d97..289daf5eb 100644
--- a/rocketpy/plots/flight_plots.py
+++ b/rocketpy/plots/flight_plots.py
@@ -400,9 +400,11 @@ def rail_buttons_forces(self):
             )
             ax1.set_xlim(
                 0,
-                self.flight.out_of_rail_time
-                if self.flight.out_of_rail_time > 0
-                else self.flight.tFinal,
+                (
+                    self.flight.out_of_rail_time
+                    if self.flight.out_of_rail_time > 0
+                    else self.flight.tFinal
+                ),
             )
             ax1.legend()
             ax1.grid(True)
@@ -431,9 +433,11 @@ def rail_buttons_forces(self):
             )
             ax2.set_xlim(
                 0,
-                self.flight.out_of_rail_time
-                if self.flight.out_of_rail_time > 0
-                else self.flight.tFinal,
+                (
+                    self.flight.out_of_rail_time
+                    if self.flight.out_of_rail_time > 0
+                    else self.flight.tFinal
+                ),
             )
             ax2.legend()
             ax2.grid(True)
@@ -557,9 +561,11 @@ def energy_data(self):
         )
         ax1.set_xlim(
             0,
-            self.flight.apogee_time
-            if self.flight.apogee_time != 0.0
-            else self.flight.t_final,
+            (
+                self.flight.apogee_time
+                if self.flight.apogee_time != 0.0
+                else self.flight.t_final
+            ),
         )
         ax1.ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
         ax1.set_title("Kinetic Energy Components")
@@ -587,9 +593,11 @@ def energy_data(self):
         )
         ax2.set_xlim(
             0,
-            self.flight.apogee_time
-            if self.flight.apogee_time != 0.0
-            else self.flight.t_final,
+            (
+                self.flight.apogee_time
+                if self.flight.apogee_time != 0.0
+                else self.flight.t_final
+            ),
         )
         ax2.ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
         ax2.set_title("Total Mechanical Energy Components")
@@ -620,9 +628,11 @@ def energy_data(self):
         )
         ax4.set_xlim(
             0,
-            self.flight.apogee_time
-            if self.flight.apogee_time != 0.0
-            else self.flight.t_final,
+            (
+                self.flight.apogee_time
+                if self.flight.apogee_time != 0.0
+                else self.flight.t_final
+            ),
         )
         ax3.ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
         ax4.set_title("Drag Absolute Power")
diff --git a/rocketpy/rocket/aero_surface.py b/rocketpy/rocket/aero_surface.py
index cea919141..a72fb74cc 100644
--- a/rocketpy/rocket/aero_surface.py
+++ b/rocketpy/rocket/aero_surface.py
@@ -1180,10 +1180,7 @@ def evaluate_geometrical_parameters(self):
             * (self.root_chord + 2 * self.tip_chord)
             * self.rocket_radius
             * self.span**2
-            + 6
-            * (self.root_chord + self.tip_chord)
-            * self.span
-            * self.rocket_radius**2
+            + 6 * (self.root_chord + self.tip_chord) * self.span * self.rocket_radius**2
         ) / 12
         roll_damping_interference_factor = 1 + (
             ((tau - λ) / (tau)) - ((1 - λ) / (tau - 1)) * np.log(tau)
@@ -1506,8 +1503,7 @@ def evaluate_geometrical_parameters(self):
                     * self.rocket_radius**2
                     * np.sqrt(-self.span**2 + self.rocket_radius**2)
                     * np.arctan(
-                        (self.span)
-                        / (np.sqrt(-self.span**2 + self.rocket_radius**2))
+                        (self.span) / (np.sqrt(-self.span**2 + self.rocket_radius**2))
                     )
                     - np.pi
                     * self.rocket_radius**2
diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index 39b0c5c7e..d72370140 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -683,9 +683,9 @@ def __init__(
                 node.time_bound = phase.TimeNodes[node_index + 1].t
                 phase.solver.t_bound = node.time_bound
                 phase.solver._lsoda_solver._integrator.rwork[0] = phase.solver.t_bound
-                phase.solver._lsoda_solver._integrator.call_args[
-                    4
-                ] = phase.solver._lsoda_solver._integrator.rwork
+                phase.solver._lsoda_solver._integrator.call_args[4] = (
+                    phase.solver._lsoda_solver._integrator.rwork
+                )
                 phase.solver.status = "running"
 
                 # Feed required parachute and discrete controller triggers
@@ -1275,9 +1275,7 @@ def udot_rail1(self, t, u, post_processing=False):
         R3 = -0.5 * rho * (free_stream_speed**2) * self.rocket.area * (drag_coeff)
 
         # Calculate Linear acceleration
-        a3 = (R3 + thrust) / M - (
-            e0**2 - e1**2 - e2**2 + e3**2
-        ) * self.env.gravity(z)
+        a3 = (R3 + thrust) / M - (e0**2 - e1**2 - e2**2 + e3**2) * self.env.gravity(z)
         if a3 > 0:
             ax = 2 * (e1 * e3 + e0 * e2) * a3
             ay = 2 * (e2 * e3 - e0 * e1) * a3
@@ -1469,9 +1467,7 @@ def u_dot(self, t, u, post_processing=False):
                         0.5 * rho * (comp_stream_speed**2) * reference_area * c_lift
                     )
                     # component lift force components
-                    lift_dir_norm = (
-                        comp_stream_vx_b**2 + comp_stream_vy_b**2
-                    ) ** 0.5
+                    lift_dir_norm = (comp_stream_vx_b**2 + comp_stream_vy_b**2) ** 0.5
                     comp_lift_xb = comp_lift * (comp_stream_vx_b / lift_dir_norm)
                     comp_lift_yb = comp_lift * (comp_stream_vy_b / lift_dir_norm)
                     # add to total lift force
@@ -1750,9 +1746,7 @@ def u_dot_generalized(self, t, u, post_processing=False):
                         0.5 * rho * (comp_stream_speed**2) * reference_area * c_lift
                     )
                     # Component lift force components
-                    lift_dir_norm = (
-                        comp_stream_vx_b**2 + comp_stream_vy_b**2
-                    ) ** 0.5
+                    lift_dir_norm = (comp_stream_vx_b**2 + comp_stream_vy_b**2) ** 0.5
                     comp_lift_xb = comp_lift * (comp_stream_vx_b / lift_dir_norm)
                     comp_lift_yb = comp_lift * (comp_stream_vy_b / lift_dir_norm)
                     # Add to total lift force
@@ -1896,8 +1890,7 @@ def u_dot_parachute(self, t, u, post_processing=False):
         freestream_z = vz
         # Determine drag force
         pseudoD = (
-            -0.5 * rho * cd_s * free_stream_speed
-            - ka * rho * 4 * np.pi * (R**2) * Rdot
+            -0.5 * rho * cd_s * free_stream_speed - ka * rho * 4 * np.pi * (R**2) * Rdot
         )
         Dx = pseudoD * freestream_x
         Dy = pseudoD * freestream_y
@@ -2486,9 +2479,7 @@ def translational_energy(self):
         # Redefine total_mass time grid to allow for efficient Function algebra
         total_mass = deepcopy(self.rocket.total_mass)
         total_mass.set_discrete_based_on_model(self.vz)
-        translational_energy = (
-            0.5 * total_mass * (self.vx**2 + self.vy**2 + self.vz**2)
-        )
+        translational_energy = 0.5 * total_mass * (self.vx**2 + self.vy**2 + self.vz**2)
         return translational_energy
 
     @funcify_method("Time (s)", "Kinetic Energy (J)", "spline", "zero")
@@ -3400,19 +3391,23 @@ def add(self, flight_phase, index=None):
                 )
                 if flight_phase.t == previous_phase.t
                 else (
-                    "Trying to add flight phase starting *together* with the one *proceeding* it. ",
-                    "This may be caused by multiple parachutes being triggered simultaneously.",
-                )
-                if flight_phase.t == next_phase.t
-                else (
-                    "Trying to add flight phase starting *before* the one *preceding* it. ",
-                    "This may be caused by multiple parachutes being triggered simultaneously",
-                    " or by having a negative parachute lag.",
-                )
-                if flight_phase.t < previous_phase.t
-                else (
-                    "Trying to add flight phase starting *after* the one *proceeding* it.",
-                    "This may be caused by multiple parachutes being triggered simultaneously.",
+                    (
+                        "Trying to add flight phase starting *together* with the one *proceeding* it. ",
+                        "This may be caused by multiple parachutes being triggered simultaneously.",
+                    )
+                    if flight_phase.t == next_phase.t
+                    else (
+                        (
+                            "Trying to add flight phase starting *before* the one *preceding* it. ",
+                            "This may be caused by multiple parachutes being triggered simultaneously",
+                            " or by having a negative parachute lag.",
+                        )
+                        if flight_phase.t < previous_phase.t
+                        else (
+                            "Trying to add flight phase starting *after* the one *proceeding* it.",
+                            "This may be caused by multiple parachutes being triggered simultaneously.",
+                        )
+                    )
                 )
             )
             self.display_warning(*warning_msg)
@@ -3420,9 +3415,7 @@ def add(self, flight_phase, index=None):
             new_index = (
                 index - 1
                 if flight_phase.t < previous_phase.t
-                else index + 1
-                if flight_phase.t > next_phase.t
-                else index
+                else index + 1 if flight_phase.t > next_phase.t else index
             )
             flight_phase.t += adjust
             self.add(flight_phase, new_index)
diff --git a/rocketpy/simulation/flight_data_importer.py b/rocketpy/simulation/flight_data_importer.py
index 3b6d3f473..edb2390d6 100644
--- a/rocketpy/simulation/flight_data_importer.py
+++ b/rocketpy/simulation/flight_data_importer.py
@@ -1,6 +1,7 @@
 """Starts with .csv or .txt file containing the rocket's collected flight data
 and build a rocketpy.Flight object from it.
 """
+
 import warnings
 from os import listdir
 from os.path import isfile, join
diff --git a/tests/test_flight_data_importer.py b/tests/test_flight_data_importer.py
index 9ce4e34ce..2dd66f3ae 100644
--- a/tests/test_flight_data_importer.py
+++ b/tests/test_flight_data_importer.py
@@ -1,5 +1,6 @@
 """Tests the FlightDataImporter class from rocketpy.simulation module.
 """
+
 import numpy as np
 
 from rocketpy.simulation import FlightDataImporter
diff --git a/tests/test_function.py b/tests/test_function.py
index 6896b01af..2ce94f691 100644
--- a/tests/test_function.py
+++ b/tests/test_function.py
@@ -49,7 +49,7 @@ def test_func_from_csv_with_header(csv_file):
     line. It tests cases where the fields are separated by quotes and without
     quotes."""
     f = Function(csv_file)
-    assert f.__repr__() == "'Function from R1 to R1 : (Scalar) → (Scalar)'"
+    assert f.__repr__() == "'Function from R1 to R1 : (time) → (value)'"
     assert np.isclose(f(0), 100)
     assert np.isclose(f(0) + f(1), 300), "Error summing the values of the function"
 
diff --git a/tests/test_genericmotor.py b/tests/test_genericmotor.py
index d9c3ff5fd..513fca40d 100644
--- a/tests/test_genericmotor.py
+++ b/tests/test_genericmotor.py
@@ -118,9 +118,7 @@ def test_generic_motor_inertia(generic_motor):
     # Tests the inertia formulation from the propellant mass
     propellant_mass = generic_motor.propellant_mass.set_discrete(*burn_time, 50).y_array
 
-    propellant_I_11 = propellant_mass * (
-        chamber_radius**2 / 4 + chamber_height**2 / 12
-    )
+    propellant_I_11 = propellant_mass * (chamber_radius**2 / 4 + chamber_height**2 / 12)
     propellant_I_22 = propellant_I_11
     propellant_I_33 = propellant_mass * (chamber_radius**2 / 2)
 
diff --git a/tests/test_tank.py b/tests/test_tank.py
index 9a6ddb6b0..6b2f03bc9 100644
--- a/tests/test_tank.py
+++ b/tests/test_tank.py
@@ -631,11 +631,7 @@ def upper_spherical_cap_centroid(radius, height=None):
     """
     if height is None:
         height = radius
-    return (
-        0.75
-        * (height**3 - 2 * height * radius**2)
-        / (height**2 - 3 * radius**2)
-    )
+    return 0.75 * (height**3 - 2 * height * radius**2) / (height**2 - 3 * radius**2)
 
 
 def tank_centroid_function(tank_radius, tank_height, zero_height=0):
@@ -775,10 +771,7 @@ def lower_spherical_cap_inertia(radius, height=None, reference=0):
     )
     inertia_y = inertia_x
     inertia_z = lower_spherical_cap_volume(radius, height) * (
-        np.pi
-        * height**3
-        * (3 * height**2 - 15 * height * radius + 20 * radius**2)
-        / 30
+        np.pi * height**3 * (3 * height**2 - 15 * height * radius + 20 * radius**2) / 30
     )
     return np.array([inertia_x, inertia_y, inertia_z])
 
diff --git a/tests/unit/test_rocket.py b/tests/unit/test_rocket.py
index 1d8cfac21..7d942ee60 100644
--- a/tests/unit/test_rocket.py
+++ b/tests/unit/test_rocket.py
@@ -207,8 +207,7 @@ def test_add_fins_assert_cp_cm_plus_fins(calisto, dimensionless_calisto, m):
         1
         + np.sqrt(
             1
-            + (2 * np.sqrt((0.12 / 2 - 0.04 / 2) ** 2 + 0.1**2) / (0.120 + 0.040))
-            ** 2
+            + (2 * np.sqrt((0.12 / 2 - 0.04 / 2) ** 2 + 0.1**2) / (0.120 + 0.040)) ** 2
         )
     )
     clalpha *= 1 + calisto.radius / (0.1 + calisto.radius)
diff --git a/tests/unit/test_solidmotor.py b/tests/unit/test_solidmotor.py
index b23dfc31e..6ae0f68ca 100644
--- a/tests/unit/test_solidmotor.py
+++ b/tests/unit/test_solidmotor.py
@@ -70,9 +70,7 @@ def test_evaluate_inertia_33_asserts_extreme_values(cesaroni_m1670):
     grain_mass = grain_vol * grain_density
 
     grain_I_33_initial = (
-        grain_mass
-        * (1 / 2.0)
-        * (grain_initial_inner_radius**2 + grain_outer_radius**2)
+        grain_mass * (1 / 2.0) * (grain_initial_inner_radius**2 + grain_outer_radius**2)
     )
 
     # not passing because I_33 is not discrete anymore