Repository: supaero-aircraft-design/FAST-GA (language: python)

Issue:
altitude_in_feet missiing from mission computation
In the computation of the descent, more precisely when computing the  calibrated airspeed for descent, an Atmosphere object is instanciated bu the altitude_in_feet attribute was not set to False which means we do not compute the descent properly


Discussion / hints:

Fix patch:
diff --git a/.codacy.yml b/.codacy.yml
index 447f84c1..ce9fab6c 100644
--- a/.codacy.yml
+++ b/.codacy.yml
@@ -3,4 +3,8 @@ exclude_paths:
   - '**/tests/**'
   - 'tests/**'
   - '**/unitary_tests/**'
-  - 'integration_tests/**'
\ No newline at end of file
+  - 'integration_tests/**'
+  - '**/conftest.py'
+  - 'docs/conf.py'
+  - '.github/**'
+  - '**.md'
\ No newline at end of file
diff --git a/src/fastga/models/aerodynamics/aero_center.py b/src/fastga/models/aerodynamics/aero_center.py
index 7d183d4e..1ac7ff27 100644
--- a/src/fastga/models/aerodynamics/aero_center.py
+++ b/src/fastga/models/aerodynamics/aero_center.py
@@ -87,10 +87,11 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         beta = math.sqrt(1.0 - mach ** 2.0)
         cl_delta_ht_cruise = cl_delta_ht / beta
 
-        # The cl_alpha_ht in the formula for the free_elevator_factor is defined with respect to the tail angle of
-        # attack, the one we compute is wth respect to the plane so it includes downwash, as a consequence we must
-        # correct it influence for this specific calculation. We will use the formula for elliptical wing as it is well
-        # known
+        # The cl_alpha_ht in the formula for the free_elevator_factor is defined with respect to
+        # the tail angle of attack, the one we compute is wth respect to the plane so it includes
+        # downwash, as a consequence we must correct it influence for this specific calculation.
+        # We will use the formula for elliptical wing as it is well known
+
         downwash_effect = 1.0 - 2.0 * cl_alpha_wing / (math.pi * aspect_ratio)
         cl_alpha_ht_ht = cl_alpha_ht / downwash_effect
         free_elevator_factor = 1.0 - (cl_delta_ht_cruise / cl_alpha_ht_ht) * (
diff --git a/src/fastga/models/aerodynamics/components/figure_digitization.py b/src/fastga/models/aerodynamics/components/figure_digitization.py
index 8edcd5ab..dc8d4f20 100644
--- a/src/fastga/models/aerodynamics/components/figure_digitization.py
+++ b/src/fastga/models/aerodynamics/components/figure_digitization.py
@@ -38,6 +38,9 @@
 KB_FLAPS = "kb_flaps.csv"
 A_DELTA_AIRFOIL = "a_delta_airfoil.csv"
 K_A_DELTA = "k_a_delta.csv"
+K_P_FLAPS = "k_p.csv"
+DELTA_CM_DELTA_CL_REF = "delta_cm_delta_cl_ref.csv"
+K_DELTA = "k_delta.csv"
 K_AR_FUSELAGE = "k_ar_fuselage.csv"
 K_VH = "k_vh.csv"
 K_CH_ALPHA = "k_ch_alpha.csv"
@@ -801,6 +804,207 @@ def k_a_delta(a_delta_airfoil, aspect_ratio) -> float:
 
         return k_a_delta
 
+    @staticmethod
+    def k_p_flaps(taper_ratio, eta_in, eta_out) -> float:
+        """
+        Roskam data to account for the partial span flaps factor on the pitch moment coefficient
+
+        :param taper_ratio: lifting surface taper ratio
+        :param eta_in: start of the control surface, in percent of the lifting surface span
+        :param eta_out: end of the control surface, in percent of the lifting surface span
+        :return k_p: partial span factor.
+        """
+
+        file = pth.join(resources.__path__[0], K_P_FLAPS)
+        db = read_csv(file)
+
+        eta_in_1_0 = FigureDigitization.interpolate_database(
+            db, "taper_1_0_X", "taper_1_0_Y", eta_in
+        )
+        eta_out_1_0 = FigureDigitization.interpolate_database(
+            db, "taper_1_0_X", "taper_1_0_Y", eta_out
+        )
+
+        eta_in_0_5 = FigureDigitization.interpolate_database(
+            db, "taper_0_5_X", "taper_0_5_Y", eta_in
+        )
+        eta_out_0_5 = FigureDigitization.interpolate_database(
+            db, "taper_0_5_X", "taper_0_5_Y", eta_out
+        )
+
+        eta_in_0_333 = FigureDigitization.interpolate_database(
+            db, "taper_0_333_X", "taper_0_333_Y", eta_in
+        )
+        eta_out_0_333 = FigureDigitization.interpolate_database(
+            db, "taper_0_333_X", "taper_0_333_Y", eta_out
+        )
+
+        eta_in_0_25 = FigureDigitization.interpolate_database(
+            db, "taper_0_25_X", "taper_0_25_Y", eta_in
+        )
+        eta_out_0_25 = FigureDigitization.interpolate_database(
+            db, "taper_0_25_X", "taper_0_25_Y", eta_out
+        )
+
+        taper_array = [0.25, 0.333, 0.5, 1.0]
+        eta_in_array = [eta_in_0_25, eta_in_0_333, eta_in_0_5, eta_in_1_0]
+        eta_out_array = [eta_out_0_25, eta_out_0_333, eta_out_0_5, eta_out_1_0]
+
+        if taper_ratio != np.clip(taper_ratio, 0.25, 1.0):
+            _LOGGER.warning("Taper ratio outside of the range in Roskam's book, value clipped")
+
+        k_p = interpolate.interp1d(taper_array, eta_out_array)(
+            np.clip(taper_ratio, 0.25, 1.0)
+        ) - interpolate.interp1d(taper_array, eta_in_array)(np.clip(taper_ratio, 0.25, 1.0))
+
+        return k_p
+
+    @staticmethod
+    def pitch_to_reference_lift(thickness_ratio: float, chord_ratio: float) -> float:
+        """
+        Roskam data to account for the ratio between the pitch moment coefficient and the
+        reference lift coefficient increment
+
+        :param thickness_ratio: thickness to chord ratio of the lifting surface
+        :param chord_ratio: chord ratio of the control surface over the lifting surface
+        :return delta_cm_delta_cl_ref: ration between the pitching moment and the reference lift
+        coefficient.
+        """
+
+        if chord_ratio != np.clip(chord_ratio, 0.05, 0.4):
+            _LOGGER.warning("Chord ratio outside of the range in Roskam's book, value clipped")
+
+        file = pth.join(resources.__path__[0], DELTA_CM_DELTA_CL_REF)
+        db = read_csv(file)
+
+        x_21 = db["TOC_21_X"]
+        y_21 = db["TOC_21_Y"]
+        errors = np.logical_or(np.isnan(x_21), np.isnan(y_21))
+        x_21 = x_21[np.logical_not(errors)].tolist()
+        y_21 = y_21[np.logical_not(errors)].tolist()
+        k_21 = interpolate.interp1d(x_21, y_21)(np.clip(chord_ratio, min(x_21), max(x_21)))
+
+        x_18 = db["TOC_18_X"]
+        y_18 = db["TOC_18_Y"]
+        errors = np.logical_or(np.isnan(x_18), np.isnan(y_18))
+        x_18 = x_18[np.logical_not(errors)].tolist()
+        y_18 = y_18[np.logical_not(errors)].tolist()
+        k_18 = interpolate.interp1d(x_18, y_18)(np.clip(chord_ratio, min(x_18), max(x_18)))
+
+        x_15 = db["TOC_15_X"]
+        y_15 = db["TOC_15_Y"]
+        errors = np.logical_or(np.isnan(x_15), np.isnan(y_15))
+        x_15 = x_15[np.logical_not(errors)].tolist()
+        y_15 = y_15[np.logical_not(errors)].tolist()
+        k_15 = interpolate.interp1d(x_15, y_15)(np.clip(chord_ratio, min(x_15), max(x_15)))
+
+        x_12 = db["TOC_12_X"]
+        y_12 = db["TOC_12_Y"]
+        errors = np.logical_or(np.isnan(x_12), np.isnan(y_12))
+        x_12 = x_12[np.logical_not(errors)].tolist()
+        y_12 = y_12[np.logical_not(errors)].tolist()
+        k_12 = interpolate.interp1d(x_12, y_12)(np.clip(chord_ratio, min(x_12), max(x_12)))
+
+        x_09 = db["TOC_09_X"]
+        y_09 = db["TOC_09_Y"]
+        errors = np.logical_or(np.isnan(x_09), np.isnan(y_09))
+        x_09 = x_09[np.logical_not(errors)].tolist()
+        y_09 = y_09[np.logical_not(errors)].tolist()
+        k_09 = interpolate.interp1d(x_09, y_09)(np.clip(chord_ratio, min(x_09), max(x_09)))
+
+        x_06 = db["TOC_06_X"]
+        y_06 = db["TOC_06_Y"]
+        errors = np.logical_or(np.isnan(x_06), np.isnan(y_06))
+        x_06 = x_06[np.logical_not(errors)].tolist()
+        y_06 = y_06[np.logical_not(errors)].tolist()
+        k_06 = interpolate.interp1d(x_06, y_06)(np.clip(chord_ratio, min(x_06), max(x_06)))
+
+        x_03 = db["TOC_03_X"]
+        y_03 = db["TOC_03_Y"]
+        errors = np.logical_or(np.isnan(x_03), np.isnan(y_03))
+        x_03 = x_03[np.logical_not(errors)].tolist()
+        y_03 = y_03[np.logical_not(errors)].tolist()
+        k_03 = interpolate.interp1d(x_03, y_03)(np.clip(chord_ratio, min(x_03), max(x_03)))
+
+        toc_array = [0.03, 0.06, 0.09, 0.12, 0.15, 0.18, 0.21]
+        k_array = [k_03, k_06, k_09, k_12, k_15, k_18, k_21]
+
+        if thickness_ratio != np.clip(thickness_ratio, 0.03, 0.4):
+            _LOGGER.warning(
+                "Thickness to chord ratio outside of the range in Roskam's book, " "value clipped"
+            )
+
+        k = interpolate.interp1d(toc_array, k_array)(np.clip(thickness_ratio, 0.03, 0.4))
+
+        return k
+
+    @staticmethod
+    def k_delta_flaps(taper_ratio: float, eta_in: float, eta_out: float) -> float:
+        """
+        Roskam data to estimate the conversion factor which accounts for partial span flaps on a
+        swept wing (c_prime/c = 1.0)
+
+        :param taper_ratio: lifting surface taper ratio
+        :param eta_in: start of the control surface, in percent of the lifting surface span
+        :param eta_out: end of the control surface, in percent of the lifting surface span
+        :return delta_k: partial span factor.
+        """
+
+        file = pth.join(resources.__path__[0], K_DELTA)
+        db = read_csv(file)
+
+        x_10 = db["X_1_0"]
+        y_10 = db["Y_1_0"]
+        errors = np.logical_or(np.isnan(x_10), np.isnan(y_10))
+        x_10 = x_10[np.logical_not(errors)].tolist()
+        y_10 = y_10[np.logical_not(errors)].tolist()
+        eta_in_1_0 = interpolate.interp1d(x_10, y_10)(np.clip(eta_in, min(x_10), max(x_10)))
+        eta_out_1_0 = interpolate.interp1d(x_10, y_10)(np.clip(eta_out, min(x_10), max(x_10)))
+
+        x_05 = db["X_0_5"]
+        y_05 = db["Y_0_5"]
+        errors = np.logical_or(np.isnan(x_05), np.isnan(y_05))
+        x_05 = x_05[np.logical_not(errors)].tolist()
+        y_05 = y_05[np.logical_not(errors)].tolist()
+        eta_in_0_5 = interpolate.interp1d(x_05, y_05)(np.clip(eta_in, min(x_05), max(x_05)))
+        eta_out_0_5 = interpolate.interp1d(x_05, y_05)(np.clip(eta_out, min(x_05), max(x_05)))
+
+        x_0333 = db["X_0_333"]
+        y_0333 = db["Y_0_333"]
+        errors = np.logical_or(np.isnan(x_0333), np.isnan(y_0333))
+        x_0333 = x_0333[np.logical_not(errors)].tolist()
+        y_0333 = y_0333[np.logical_not(errors)].tolist()
+        eta_in_0_333 = interpolate.interp1d(x_0333, y_0333)(
+            np.clip(eta_in, min(x_0333), max(x_0333))
+        )
+        eta_out_0_333 = interpolate.interp1d(x_0333, y_0333)(
+            np.clip(eta_out, min(x_0333), max(x_0333))
+        )
+
+        x_02 = db["X_0_2"]
+        y_02 = db["Y_0_2"]
+        errors = np.logical_or(np.isnan(x_02), np.isnan(y_02))
+        x_02 = x_02[np.logical_not(errors)].tolist()
+        y_02 = y_02[np.logical_not(errors)].tolist()
+        eta_in_0_2 = interpolate.interp1d(x_02, y_02)(np.clip(eta_in, min(x_02), max(x_02)))
+        eta_out_0_2 = interpolate.interp1d(x_02, y_02)(np.clip(eta_out, min(x_02), max(x_02)))
+
+        taper_array = [0.2, 0.333, 0.5, 1.0]
+        eta_in_array = [eta_in_0_2, eta_in_0_333, eta_in_0_5, eta_in_1_0]
+        eta_out_array = [eta_out_0_2, eta_out_0_333, eta_out_0_5, eta_out_1_0]
+
+        if taper_ratio != np.clip(taper_ratio, 0.2, 1.0):
+            _LOGGER.warning("Taper ratio outside of the range in Roskam's book, value clipped")
+
+        k_delta_in = interpolate.interp1d(taper_array, eta_in_array)(np.clip(taper_ratio, 0.2, 1.0))
+        k_delta_out = interpolate.interp1d(taper_array, eta_out_array)(
+            np.clip(taper_ratio, 0.2, 1.0)
+        )
+
+        k_delta = k_delta_out - k_delta_in
+
+        return k_delta
+
     @staticmethod
     def k_ar_fuselage(taper_ratio, span, avg_fuselage_depth) -> float:
         """
@@ -812,7 +1016,7 @@ def k_ar_fuselage(taper_ratio, span, avg_fuselage_depth) -> float:
         :param avg_fuselage_depth: average fuselage depth (diameter if fuselage considered
         circular), in m
         :return k_ar_fuselage: correction factor to account for the end plate effect of the fuselage
-         on effective VTP AR
+         on effective VTP AR.
         """
 
         file = pth.join(resources.__path__[0], K_AR_FUSELAGE)
@@ -1133,3 +1337,18 @@ def k_fus(root_quarter_chord_position_ratio) -> float:
         )
 
         return k_fus
+
+    @staticmethod
+    def interpolate_database(database, tag_x: str, tag_y: str, input_x: float):
+
+        database_x = database[tag_x]
+        database_y = database[tag_y]
+        errors = np.logical_or(np.isnan(database_x), np.isnan(database_x))
+        database_x = database_x[np.logical_not(errors)].tolist()
+        database_y = database_y[np.logical_not(errors)].tolist()
+
+        output_y = interpolate.interp1d(database_x, database_y)(
+            np.clip(input_x, min(database_x), max(database_x))
+        )
+
+        return output_y
diff --git a/src/fastga/models/aerodynamics/components/high_lift_aero.py b/src/fastga/models/aerodynamics/components/high_lift_aero.py
index 90144064..e86e89a1 100644
--- a/src/fastga/models/aerodynamics/components/high_lift_aero.py
+++ b/src/fastga/models/aerodynamics/components/high_lift_aero.py
@@ -219,15 +219,48 @@ def _get_flaps_delta_cm(self, inputs, flap_angle: float, mach: float) -> float:
         :return: increment of moment coefficient
         """
 
-        wing_taper_ratio = float(inputs["data:geometry:wing:taper_ratio"])
+        taper_ratio_wing = inputs["data:geometry:wing:taper_ratio"]
+        cl_alpha_airfoil_wing = inputs["data:aerodynamics:wing:airfoil:CL_alpha"]
+        span_wing = inputs["data:geometry:wing:span"]
+        y1_wing = inputs["data:geometry:fuselage:maximum_width"] / 2.0
+        y2_wing = inputs["data:geometry:wing:root:y"]
+        flap_span_ratio = inputs["data:geometry:flap:span_ratio"]
+        taper_ratio_wing = float(inputs["data:geometry:wing:taper_ratio"])
+        aspect_ratio_wing = float(inputs["data:geometry:wing:aspect_ratio"])
+        flap_chord_ratio = float(inputs["data:geometry:flap:chord_ratio"])
+        wing_thickness_ratio = float(inputs["data:geometry:wing:thickness_ratio"])
+        sweep_25 = float(inputs["data:geometry:wing:sweep_25"]) * np.pi / 180.0
+
+        beta_ref = np.sqrt(1.0 - mach ** 2.0)
+        k = cl_alpha_airfoil_wing / (2.0 * np.pi)
+        cl_alpha_ref = (
+            2.0 * np.pi * 6.0 / (2.0 + np.sqrt((36.0 * beta_ref ** 2.0 / k ** 2.0 + 4.0)))
+        )
+
+        # First we need to compute the increment in the lift coefficient for the reference case
+        delta_cl_2d_ref = self._compute_delta_cl_airfoil_2d(inputs, flap_angle, mach)
+        eta_in_ref = 0.0
+        eta_out_ref = 1.0
+        kb_ref = self.k_b_flaps(eta_in_ref, eta_out_ref, taper_ratio_wing)
+        a_delta_flap_ref = self.a_delta_airfoil(float(inputs["data:geometry:flap:chord_ratio"]))
+        k_a_delta_ref = self.k_a_delta(a_delta_flap_ref, 6.0)
+        delta_cl_ref = (
+            kb_ref * delta_cl_2d_ref * (cl_alpha_ref / cl_alpha_airfoil_wing) * k_a_delta_ref
+        )
+
+        # Now we can compute the coefficient
+        eta_in = float(y1_wing / (span_wing / 2.0))
+        eta_out = float(
+            ((y2_wing - y1_wing) + flap_span_ratio * (span_wing / 2.0 - y2_wing))
+            / (span_wing / 2.0 - y2_wing)
+        )
+        k_delta = self.k_delta_flaps(taper_ratio_wing, eta_in, eta_out)
+        k_p = self.k_p_flaps(taper_ratio_wing, eta_in, eta_out)
+        delta_cm_delta_cl_ref = self.pitch_to_reference_lift(wing_thickness_ratio, flap_chord_ratio)
 
-        # Method from Roskam (sweep=0, flaps 60%, simple slotted and not extensible, at 25% MAC, cf/c+0.25)
-        k_p = interpolate.interp1d([0.0, 0.2, 0.33, 0.5, 1.0], [0.65, 0.75, 0.7, 0.63, 0.5])
-        # k_p: Figure 8.105, interpolated function of taper ratio (span ratio fixed)
-        delta_cl_flap = self._get_flaps_delta_cl(inputs, flap_angle, mach)[0]
         delta_cm_flap = (
-            k_p(min(max(0.0, wing_taper_ratio), 1.0)) * (-0.27) * delta_cl_flap
-        )  # -0.27: Figure 8.106
+            k_delta * aspect_ratio_wing / 1.5 * np.tan(sweep_25) + k_p * delta_cm_delta_cl_ref
+        ) * delta_cl_ref
 
         return delta_cm_flap
 
diff --git a/src/fastga/models/aerodynamics/components/resources/delta_cm_delta_cl_ref.csv b/src/fastga/models/aerodynamics/components/resources/delta_cm_delta_cl_ref.csv
new file mode 100644
index 00000000..2e5c48c9
--- /dev/null
+++ b/src/fastga/models/aerodynamics/components/resources/delta_cm_delta_cl_ref.csv
@@ -0,0 +1,16 @@
+TOC_21_X,TOC_21_Y,TOC_18_X,TOC_18_Y,TOC_15_X,TOC_15_Y,TOC_12_X,TOC_12_Y,TOC_09_X,TOC_09_Y,TOC_06_X,TOC_06_Y,TOC_03_X,TOC_03_Y
+0.04947589098532494,-0.3010738255033557,0.048951781970649894,-0.32120805369127514,0.05,-0.3542281879194631,0.04947589098532494,-0.38966442953020136,0.047903563941299784,-0.430738255033557,0.047903563941299784,-0.4677852348993289,0.04842767295597483,-0.5076510067114094
+0.0741090146750524,-0.2934228187919463,0.07515723270440251,-0.31154362416107384,0.07358490566037736,-0.3409395973154362,0.07463312368972745,-0.37395973154362416,0.0730607966457023,-0.4093959731543624,0.07253668763102725,-0.44281879194630874,0.0741090146750524,-0.4774496644295302
+0.09926624737945491,-0.2873825503355705,0.09874213836477988,-0.3030872483221476,0.09926624737945491,-0.329261744966443,0.09926624737945491,-0.3578523489932886,0.09979035639412998,-0.38845637583892617,0.09874213836477988,-0.41825503355704696,0.09926624737945491,-0.4468456375838926
+0.12494758909853249,-0.28013422818791944,0.12389937106918239,-0.2934228187919463,0.12442348008385742,-0.3175838926174497,0.12442348008385742,-0.3433557046979866,0.12442348008385742,-0.3687248322147651,0.12442348008385742,-0.39248322147651005,0.12389937106918239,-0.4174496644295302
+0.150104821802935,-0.2724832214765101,0.14958071278825996,-0.2857718120805369,0.15062893081761006,-0.30590604026845636,0.150104821802935,-0.327248322147651,0.150104821802935,-0.34818791946308725,0.14958071278825996,-0.3691275167785235,0.14905660377358487,-0.3912751677852349
+0.1752620545073375,-0.2664429530201342,0.17473794549266244,-0.27691275167785234,0.17473794549266244,-0.29422818791946315,0.1752620545073375,-0.31275167785234903,0.1742138364779874,-0.3300671140939597,0.17473794549266244,-0.34697986577181206,0.1752620545073375,-0.3642953020134228
+0.19937106918238992,-0.26040268456375837,0.19989517819706498,-0.2696644295302013,0.19989517819706498,-0.28456375838926173,0.19989517819706498,-0.2978523489932886,0.19989517819706498,-0.31234899328859056,0.19989517819706498,-0.3268456375838926,0.19937106918238992,-0.33932885906040267
+0.22557651991614253,-0.2543624161073826,0.2250524109014675,-0.262013422818792,0.2245283018867924,-0.2732885906040269,0.2250524109014675,-0.2857718120805369,0.2250524109014675,-0.2954362416107383,0.22400419287211737,-0.3079194630872483,0.2245283018867924,-0.3175838926174497
+0.2496855345911949,-0.24832214765100669,0.25020964360586995,-0.25516778523489936,0.24916142557651988,-0.2640268456375839,0.25020964360586995,-0.27288590604026847,0.2496855345911949,-0.28013422818791944,0.2496855345911949,-0.2893959731543624,0.25073375262054504,-0.2986577181208054
+0.2748427672955974,-0.24389261744966445,0.2748427672955974,-0.2487248322147651,0.2753668763102725,-0.25476510067114094,0.2753668763102725,-0.262013422818792,0.2743186582809224,-0.2684563758389262,0.2748427672955974,-0.27489932885906043,0.2753668763102725,-0.28255033557046977
+0.2994758909853249,-0.23865771812080533,0.29895178197064987,-0.24268456375838926,0.29895178197064987,-0.24671140939597314,0.3,-0.25194630872483226,0.3,-0.2567785234899329,0.3,-0.262013422818792,0.29895178197064987,-0.2684563758389262
+0.32463312368972747,-0.23382550335570468,0.3251572327044025,-0.23664429530201342,0.32463312368972747,-0.23946308724832216,0.32463312368972747,-0.24268456375838926,0.32463312368972747,-0.24671140939597314,0.32463312368972747,-0.2503355704697987,0.32463312368972747,-0.25516778523489936
+0.3497903563941299,-0.2297986577181208,0.34926624737945483,-0.23140939597315435,0.3497903563941299,-0.23261744966442954,0.34926624737945483,-0.23543624161073823,0.3497903563941299,-0.23825503355704697,0.3497903563941299,-0.24026845637583893,0.3497903563941299,-0.24308724832214768
+0.37547169811320746,-0.22577181208053693,0.37494758909853243,-0.22778523489932884,0.37547169811320746,-0.22778523489932884,0.37547169811320746,-0.22899328859060403,0.37494758909853243,-0.2297986577181208,0.37494758909853243,-0.23140939597315435,0.37494758909853243,-0.23261744966442954
+0.40010482180293494,-0.22255033557046983,0.40062893081761003,-0.2221476510067114,0.40010482180293494,-0.22255033557046983,0.3995807127882599,-0.2233557046979866,0.40062893081761003,-0.2233557046979866,0.40010482180293494,-0.22255033557046983,0.40010482180293494,-0.22255033557046983
diff --git a/src/fastga/models/aerodynamics/components/resources/k_delta.csv b/src/fastga/models/aerodynamics/components/resources/k_delta.csv
new file mode 100644
index 00000000..09b6175a
--- /dev/null
+++ b/src/fastga/models/aerodynamics/components/resources/k_delta.csv
@@ -0,0 +1,22 @@
+X_1_0,Y_1_0,X_0_5,Y_0_5,X_0_333,Y_0_333,X_0_2,Y_0_2
+0,0,0,0,0,0,0,0
+0.10099889,0.023051771,0.002599388,0.000359897,0.032428356,0.009922078,0.02967033,0.01015873
+0.136514983,0.029918256,0.036850153,0.010025707,0.076772247,0.019688312,0.064835165,0.02
+0.203107658,0.041689373,0.101070336,0.024421594,0.102111614,0.025818182,0.103296703,0.027195767
+0.266370699,0.049754768,0.12675841,0.029974293,0.123227753,0.029974026,0.117582418,0.02994709
+0.305216426,0.053896458,0.181345566,0.039742931,0.178129713,0.03974026,0.171428571,0.03957672
+0.377358491,0.059673025,0.20382263,0.043239075,0.204524887,0.043688312,0.202197802,0.044232804
+0.405105438,0.06119891,0.303363914,0.055167095,0.253092006,0.049818182,0.252747253,0.04994709
+0.448390677,0.062724796,0.347247706,0.057943445,0.304826546,0.054701299,0.304395604,0.053756614
+0.48390677,0.062615804,0.401834862,0.060205656,0.346003017,0.057194805,0.357142857,0.055767196
+0.503884573,0.062288828,0.441437309,0.060616967,0.405128205,0.058441558,0.405494505,0.056190476
+0.56936737,0.06,0.50351682,0.059691517,0.464253394,0.057506494,0.504395604,0.052910053
+0.602663707,0.058038147,0.550611621,0.057532134,0.505429864,0.05625974,0.551648352,0.04962963
+0.658157603,0.054332425,0.603058104,0.054138817,0.564555053,0.053142857,0.604395604,0.044867725
+0.703662597,0.049754768,0.649082569,0.050437018,0.601508296,0.05012987,0.658241758,0.039470899
+0.786903441,0.039836512,0.700458716,0.044987147,0.704977376,0.040779221,0.701098901,0.035026455
+0.802441731,0.038092643,0.8,0.032956298,0.803167421,0.028311688,0.748351648,0.02952381
+0.852386238,0.029482289,0.890978593,0.019280206,0.866515837,0.019168831,0.8,0.023068783
+0.901220866,0.020871935,0.900611621,0.017532134,0.899245852,0.014597403,0.827472527,0.01978836
+0.958934517,0.009318801,0.9999999,0.00001,0.995324284,0.0001,0.897802198,0.011428571
+1,0,1,0,1,0,1,0
diff --git a/src/fastga/models/aerodynamics/components/resources/k_p.csv b/src/fastga/models/aerodynamics/components/resources/k_p.csv
new file mode 100644
index 00000000..458659e1
--- /dev/null
+++ b/src/fastga/models/aerodynamics/components/resources/k_p.csv
@@ -0,0 +1,12 @@
+taper_1_0_X,taper_1_0_Y,taper_0_5_X,taper_0_5_Y,taper_0_333_X,taper_0_333_Y,taper_0_25_X,taper_0_25_Y
+0,0,0,0,0,0,0,0
+0.10222222222222227,0.09739130434782628,0.10222222222222227,0.1617391304347826,0.10044444444444448,0.20347826086956533,0.10222222222222227,0.231304347826087
+0.2,0.1982608695652175,0.19822222222222224,0.3165217391304349,0.19822222222222224,0.3930434782608696,0.19822222222222224,0.44000000000000006
+0.2995555555555555,0.3008695652173914,0.2995555555555555,0.462608695652174,0.2995555555555555,0.5565217391304348,0.2977777777777778,0.6173913043478262
+0.39733333333333326,0.4017391304347826,0.3991111111111111,0.5930434782608696,0.3955555555555555,0.6973913043478261,0.39733333333333326,0.7617391304347827
+0.4968888888888889,0.5060869565217391,0.4968888888888889,0.7026086956521739,0.4968888888888889,0.8104347826086957,0.4968888888888889,0.8817391304347827
+0.5964444444444443,0.6034782608695652,0.5946666666666665,0.7965217391304348,0.5946666666666665,0.9008695652173914,0.5946666666666665,0.9686956521739131
+0.6924444444444444,0.7060869565217391,0.6942222222222221,0.871304347826087,0.6942222222222221,0.9704347826086956,0.6942222222222221,1.0330434782608697
+0.7919999999999998,0.8052173913043479,0.7937777777777777,0.942608695652174,0.7919999999999998,1.022608695652174,0.7937777777777777,1.08
+0.8933333333333333,0.9043478260869565,0.8915555555555554,1.0017391304347827,0.8897777777777776,1.066086956521739,0.8897777777777776,1.1113043478260871
+0.9875555555555553,1.0034782608695654,0.9875555555555553,1.0556521739130436,0.9857777777777776,1.1008695652173914,0.9875555555555553,1.128695652173913
diff --git a/src/fastga/models/aerodynamics/external/openvsp/compute_aero.py b/src/fastga/models/aerodynamics/external/openvsp/compute_aero.py
index eb36e611..c1a2ec7b 100644
--- a/src/fastga/models/aerodynamics/external/openvsp/compute_aero.py
+++ b/src/fastga/models/aerodynamics/external/openvsp/compute_aero.py
@@ -136,7 +136,7 @@ def compute(self, inputs, outputs):
         _LOGGER.debug("Entering aerodynamic computation")
 
         # Check AOA input is float
-        if not (type(INPUT_AOA) == float):
+        if not isinstance(INPUT_AOA, float):
             raise TypeError("INPUT_AOA should be a float!")
 
         if self.options["low_speed_aero"]:
@@ -172,7 +172,7 @@ def compute(self, inputs, outputs):
                 "data:aerodynamics:aircraft:mach_interpolation:CL_alpha_vector"
             ] = cl_alpha_interp
 
-        # Defining outputs -----------------------------------------------------------------------------
+        # Defining outputs
         if self.options["low_speed_aero"]:
             outputs["data:aerodynamics:wing:low_speed:CL0_clean"] = cl_0_wing
             outputs["data:aerodynamics:wing:low_speed:CL_alpha"] = cl_alpha_wing
diff --git a/src/fastga/models/aerodynamics/external/openvsp/compute_aero_slipstream.py b/src/fastga/models/aerodynamics/external/openvsp/compute_aero_slipstream.py
index d4a2dd77..12c254b6 100644
--- a/src/fastga/models/aerodynamics/external/openvsp/compute_aero_slipstream.py
+++ b/src/fastga/models/aerodynamics/external/openvsp/compute_aero_slipstream.py
@@ -149,8 +149,9 @@ def compute(self, inputs, outputs):
         atm = Atmosphere(altitude, altitude_in_feet=False)
         velocity = mach * atm.speed_of_sound
 
-        # We need to compute the AOA for which the most constraining Delta_Cl due to slipstream will appear, this
-        # is taken as the angle for which the clean wing is at its max angle of attack
+        # We need to compute the AOA for which the most constraining Delta_Cl due to slipstream
+        # will appear, this is taken as the angle for which the clean wing is at its max angle of
+        # attack
 
         alpha_max = (cl_max_clean - cl0) / cl_alpha
 
diff --git a/src/fastga/models/aerodynamics/external/vlm/compute_aero.py b/src/fastga/models/aerodynamics/external/vlm/compute_aero.py
index 91b7bd6d..46b76992 100644
--- a/src/fastga/models/aerodynamics/external/vlm/compute_aero.py
+++ b/src/fastga/models/aerodynamics/external/vlm/compute_aero.py
@@ -247,7 +247,7 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         _LOGGER.debug("Entering aerodynamic computation")
 
         # Check AOA input is float
-        if not (type(INPUT_AOA) == float):
+        if not isinstance(INPUT_AOA, float):
             raise TypeError("INPUT_AOA should be a float!")
 
         if self.options["low_speed_aero"]:
@@ -282,7 +282,7 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
                     inputs, outputs, INPUT_AOA, altitude, mach
                 )
 
-        # Defining outputs -----------------------------------------------------------------------------
+        # Defining outputs
         if self.options["low_speed_aero"]:
             outputs["data:aerodynamics:wing:low_speed:CL0_clean"] = cl_0_wing
             outputs["data:aerodynamics:wing:low_speed:CL_alpha"] = cl_alpha_wing
diff --git a/src/fastga/models/aerodynamics/external/xfoil/xfoil_polar.py b/src/fastga/models/aerodynamics/external/xfoil/xfoil_polar.py
index 921c3f19..0ee9a9ea 100644
--- a/src/fastga/models/aerodynamics/external/xfoil/xfoil_polar.py
+++ b/src/fastga/models/aerodynamics/external/xfoil/xfoil_polar.py
@@ -108,8 +108,8 @@ def compute(self, inputs, outputs):
         mach = round(float(inputs["xfoil:mach"]) * 1e4) / 1e4
         reynolds = round(float(inputs["xfoil:reynolds"]))
 
-        # Search if data already stored for this profile and mach with reynolds values bounding current value.
-        # If so, use linear interpolation with the nearest upper/lower reynolds
+        # Search if data already stored for this profile and mach with reynolds values bounding
+        # current value. If so, use linear interpolation with the nearest upper/lower reynolds
         no_file = True
         data_saved = None
         interpolated_result = None
@@ -202,7 +202,8 @@ def compute(self, inputs, outputs):
                         upper_value = np.array(
                             np.matrix(upper_values.loc[label, index_upper_reynolds].to_numpy()[0])
                         ).ravel()
-                        # If values relative to alpha vector, performs interpolation with shared vector
+                        # If values relative to alpha vector, performs interpolation with shared
+                        # vector
                         if np.size(lower_value) == len(alpha_lower):
                             lower_value = np.interp(
                                 alpha_shared, np.array(alpha_lower), lower_value
@@ -219,7 +220,7 @@ def compute(self, inputs, outputs):
             if result_folder_path != "":
                 os.makedirs(result_folder_path, exist_ok=True)
 
-            # Pre-processing (populating temp directory) -----------------------------------------------
+            # Pre-processing (populating temp directory)
             # XFoil exe
             tmp_directory = self._create_tmp_directory()
             if self.options[OPTION_XFOIL_EXE_PATH]:
@@ -263,14 +264,15 @@ def compute(self, inputs, outputs):
                 ALPHA_STEP,
             )
 
-            # Run XFOIL --------------------------------------------------------------------------------
+            # Run XFOIL
             self.options["external_input_files"] = [self.stdin, tmp_profile_file_path]
             self.options["external_output_files"] = [tmp_result_file_path]
             # noinspection PyBroadException
             try:
                 super().compute(inputs, outputs)
                 result_array_p = self._read_polar(tmp_result_file_path)
-            except:  # catch the error and try to read result file for non-convergence on higher angles
+            except:
+                # catch the error and try to read result file for non-convergence on higher angles
                 e = sys.exc_info()[1]
                 try:
                     result_array_p = self._read_polar(tmp_result_file_path)
@@ -296,14 +298,15 @@ def compute(self, inputs, outputs):
                 try:
                     super().compute(inputs, outputs)
                     result_array_n = self._read_polar(tmp_result_file_path)
-                except:  # catch the error and try to read result file for non-convergence on higher angles
+                except:
+                    # catch the error and try to read result file for non-convergence on higher angles
                     e = sys.exc_info()[1]
                     try:
                         result_array_n = self._read_polar(tmp_result_file_path)
                     except:
                         raise TimeoutError("<p>Error: %s</p>" % e)
 
-            # Post-processing --------------------------------------------------------------------------
+            # Post-processing
             if self.options[OPTION_COMP_NEG_AIR_SYM]:
                 cl_max_2d, error = self._get_max_cl(result_array_p["alpha"], result_array_p["CL"])
                 # noinspection PyUnboundLocalVariable
@@ -353,7 +356,7 @@ def compute(self, inputs, outputs):
                 cm.extend(additional_zeros)
                 cm = np.array(cm)
 
-            # Save results to defined path -------------------------------------------------------------
+            # Save results to defined path
             if not error:
                 results = [
                     np.array(mach),
@@ -390,7 +393,7 @@ def compute(self, inputs, outputs):
                         "%s folder!" % local_resources.__path__[0]
                     )
 
-            # Getting output files if needed -------------------------------------------------------
+            # Getting output files if needed
             if self.options[OPTION_RESULT_FOLDER_PATH] != "":
                 if pth.exists(tmp_result_file_path):
                     polar_file_path = pth.join(
@@ -405,8 +408,8 @@ def compute(self, inputs, outputs):
                 if pth.exists(self.stderr):
                     stderr_file_path = pth.join(result_folder_path, _STDERR_FILE_NAME)
                     shutil.move(self.stderr, stderr_file_path)
-            # Try to delete the temp directory, if process not finished correctly try to close files before removing
-            # directory for second attempt
+            # Try to delete the temp directory, if process not finished correctly try to
+            # close files before removing directory for second attempt
             # noinspection PyBroadException
             try:
                 tmp_directory.cleanup()
@@ -530,7 +533,8 @@ def _get_max_cl(self, alpha: np.ndarray, lift_coeff: np.ndarray) -> Tuple[float,
 
         :param alpha:
         :param lift_coeff: CL
-        :return: max CL within +/- 0.3 around linear zone if enough alpha computed, or default value otherwise
+        :return: max CL within +/- 0.3 around linear zone if enough alpha computed, or default value
+        otherwise
         """
         alpha_range = self.options[OPTION_ALPHA_END] - self.options[OPTION_ALPHA_START]
         if len(alpha) > 2:
@@ -573,7 +577,7 @@ def _get_min_cl(self, alpha: np.ndarray, lift_coeff: np.ndarray) -> Tuple[float,
                 return min(lift_coeff[delta <= 0.3]), False
 
         _LOGGER.warning(
-            "2D CL min not found, less than 40% of angle range computed: using default value %f",
+            "2D CL min not found, less than 40%% of angle range computed: using default value %f",
             DEFAULT_2D_CL_MIN,
         )
         return DEFAULT_2D_CL_MIN, True
diff --git a/src/fastga/models/geometry/geom_components/fuselage/components/compute_fuselage_dimensions.py b/src/fastga/models/geometry/geom_components/fuselage/components/compute_fuselage_dimensions.py
index bc3aab90..26adae1a 100644
--- a/src/fastga/models/geometry/geom_components/fuselage/components/compute_fuselage_dimensions.py
+++ b/src/fastga/models/geometry/geom_components/fuselage/components/compute_fuselage_dimensions.py
@@ -56,7 +56,8 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
 class ComputeFuselageGeometryCabinSizingFD(ExplicitComponent):
     # TODO: Document equations. Cite sources
     """
-    Geometry of fuselage - Cabin is sized based on layout (seats, aisle...) and HTP/VTP position (Fixed tail Distance).
+    Geometry of fuselage - Cabin is sized based on layout (seats, aisle...) and HTP/VTP position
+    (Fixed tail Distance).
     """
 
     def __init__(self, **kwargs):
@@ -156,8 +157,9 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
             lav = propulsion_length + spinner_length
         else:
             lav = 1.40 * h_f
-            # Used to be 1.7, supposedly as an A320 according to FAST legacy. Results on the BE76 tend to say it is
-            # around 1.40, though it varies a lot depending on the airplane and its use
+            # Used to be 1.7, supposedly as an A320 according to FAST legacy. Results on the BE76
+            # tend to say it is around 1.40, though it varies a lot depending on the airplane and
+            # its use
         # Calculate fuselage length
         fus_length = fa_length + max(ht_lp + 0.75 * ht_length, vt_lp + 0.75 * vt_length)
         plane_length = fa_length + max(
@@ -181,7 +183,8 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
 class ComputeFuselageGeometryCabinSizingFL(ExplicitComponent):
     # TODO: Document equations. Cite sources
     """
-    Geometry of fuselage - Cabin is sized based on layout (seats, aisle...) and additional rear length (Fixed Length).
+    Geometry of fuselage - Cabin is sized based on layout (seats, aisle...) and additional rear
+    length (Fixed Length).
     """
 
     def __init__(self, **kwargs):
diff --git a/src/fastga/models/geometry/profiles/profile.py b/src/fastga/models/geometry/profiles/profile.py
index f8a726de..31452fef 100644
--- a/src/fastga/models/geometry/profiles/profile.py
+++ b/src/fastga/models/geometry/profiles/profile.py
@@ -28,19 +28,19 @@
 
 
 class Profile:
-    """Class for managing 2D wing profiles"""
+    """Class for managing 2D wing profiles."""
 
-    # pylint: disable=invalid-name  # X and Z are valid names in this context
+    # pylint: disable=invalid-name
+    # X and Z are valid names in this context
 
     def __init__(self, chord_length: float = 0.0):
 
         self._rel_mean_line_and_thickness = pd.DataFrame(columns=[X, Z, THICKNESS])
         """
         Data of mean line and thickness, computed after inputs of :meth:`set_points`_.
-        
         DataFrame keys are 'x', 'z' and 'thickness'.
         - 'x' and 'z' are relative to chord_length
-        - 'thickness' is relative to max thickness (and given according to 'x')
+        - 'thickness' is relative to max thickness (and given according to 'x').
         """
 
         self.chord_length: float = chord_length
diff --git a/src/fastga/models/handling_qualities/tail_sizing/compute_balked_landing_limit.py b/src/fastga/models/handling_qualities/tail_sizing/compute_balked_landing_limit.py
index 1fe6c71c..282c9e92 100644
--- a/src/fastga/models/handling_qualities/tail_sizing/compute_balked_landing_limit.py
+++ b/src/fastga/models/handling_qualities/tail_sizing/compute_balked_landing_limit.py
@@ -116,9 +116,9 @@ def found_cl_repartition(inputs, load_factor, mass, dynamic_pressure, x_cg):
         inv_a = np.linalg.inv(a)
         CL = np.dot(inv_a, b)
 
-        # Now that we have the Cl on both lifting surfaces we need to find the corresponding aircraft angle of attack
-        # and elevator deflection to see if the equilibrium is possible, but first we must remove the effect of HLD
-        # on the wing
+        # Now that we have the Cl on both lifting surfaces we need to find the corresponding
+        # aircraft angle of attack and elevator deflection to see if the equilibrium is possible,
+        # but first we must remove the effect of HLD on the wing
 
         Cl_corrected_1 = float(CL[0] - (cl_flaps + cl0_wing))
         Cl_corrected_2 = float(CL[1])
diff --git a/src/fastga/models/load_analysis/wing/aerodynamic_loads.py b/src/fastga/models/load_analysis/wing/aerodynamic_loads.py
index 4ca3262c..d3685cf9 100644
--- a/src/fastga/models/load_analysis/wing/aerodynamic_loads.py
+++ b/src/fastga/models/load_analysis/wing/aerodynamic_loads.py
@@ -1,5 +1,6 @@
 """
-Computes the aerodynamic loads on the wing of the aircraft in the most stringent case according to aerostructural loads.
+Computes the aerodynamic loads on the wing of the aircraft in the most stringent case
+according to aerostructural loads.
 """
 
 #  This file is part of FAST : A framework for rapid Overall Aircraft Design
@@ -133,8 +134,7 @@ def setup(self):
 
     def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
 
-        # STEP 1/XX - DEFINE OR CALCULATE INPUT DATA FOR LOAD COMPUTATION ##############################################
-        ################################################################################################################
+        # STEP 1/XX - DEFINE OR CALCULATE INPUT DATA FOR LOAD COMPUTATION
 
         y_vector = inputs["data:aerodynamics:wing:low_speed:Y_vector"]
         y_vector_slip = inputs["data:aerodynamics:slipstream:wing:cruise:prop_on:Y_vector"]
@@ -156,11 +156,12 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
 
         cruise_v_tas = inputs["data:TLAR:v_cruise"]
 
-        # STEP 2/XX - DELETE THE ADDITIONAL ZEROS WE HAD TO PUT TO FIT OPENMDAO AND ADD A POINT AT THE ROOT (Y=0) AND AT
-        # THE VERY TIP (Y=SPAN/2) TO GET THE WHOLE SPAN OF THE WING IN THE INTERPOLATION WE WILL DO LATER ##############
+        # STEP 2/XX - DELETE THE ADDITIONAL ZEROS WE HAD TO PUT TO FIT OPENMDAO AND ADD A POINT
+        # AT THE ROOT (Y=0) AND AT THE VERY TIP (Y=SPAN/2) TO GET THE WHOLE SPAN OF THE WING IN
+        # THE INTERPOLATION WE WILL DO LATER
 
-        # We delete the zeros we had to add to fit the size we set in the aerodynamics module and add the physic
-        # extrema that are missing, the root and the full span,
+        # We delete the zeros we had to add to fit the size we set in the aerodynamics module and
+        # add the physic extrema that are missing, the root and the full span,
         y_vector = AerostructuralLoad.delete_additional_zeros(y_vector)
         y_vector_slip = AerostructuralLoad.delete_additional_zeros(y_vector_slip)
         cl_vector = AerostructuralLoad.delete_additional_zeros(cl_vector, len(y_vector))
@@ -182,18 +183,20 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         cl_vector_slip = np.append(cl_vector_slip, 0.0)
         chord_vector = np.append(chord_vector, tip_chord)
 
-        # To get the same y_vector array as in the aerostructural computation, the only important part here is the
-        # location of the y samples so we don't need to register the structural mass array
+        # To get the same y_vector array as in the aerostructural computation, the only important
+        # part here is the location of the y samples so we don't need to register the structural
+        # mass array
         y_vector, _ = self.compute_relief_force(
             inputs, y_vector_orig, chord_vector, 0.0, 0.0, False
         )
 
-        # STEP 3/XX - WE COMPUTE THE BASELINE LIFT AND SCALE IT UP ACCORDING TO THE MOST CONSTRAINING CASE IDENTIFIED IN
-        # THE AEROSTRUCTURAL ANALYSIS ##################################################################################
+        # STEP 3/XX - WE COMPUTE THE BASELINE LIFT AND SCALE IT UP ACCORDING TO THE MOST
+        # CONSTRAINING CASE IDENTIFIED IN THE AEROSTRUCTURAL ANALYSIS
 
-        # Now we identify the constraint that gives the highest lift distribution in amplitude which is linked with
-        # the highest absolute load factor. From there we can recompute the equilibrium and get the lift distribution
-        # in the most stringent case which will be what we will plot
+        # Now we identify the constraint that gives the highest lift distribution in amplitude
+        # which is linked with the highest absolute load factor. From there we can recompute the
+        # equilibrium and get the lift distribution in the most stringent case which will be what
+        # we will plot
         if load_factor_shear > load_factor_rbm:
             mass = inputs["data:loads:max_shear:mass"]
             load_factor = load_factor_shear
@@ -213,8 +216,8 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         cl_s_slip_actual = cl_s_slip * (v_ref / cruise_v_tas) ** 2.0
         lift_distribution = (cl_s_actual + cl_s_slip_actual) * dynamic_pressure
 
-        # STEP 4/XX - WE ADD ZEROS AT THE END OF THE RESULT LIFT DISTRIBUTION TO FIT THE FORMAT IMPOSED BY OPENMDAO ####
-        ################################################################################################################
+        # STEP 4/XX - WE ADD ZEROS AT THE END OF THE RESULT LIFT DISTRIBUTION TO FIT THE FORMAT
+        # IMPOSED BY OPENMDAO
 
         additional_zeros = np.zeros(SPAN_MESH_POINT_LOADS - len(y_vector))
         lift_distribution_outputs = np.concatenate([lift_distribution, additional_zeros])
diff --git a/src/fastga/models/load_analysis/wing/aerostructural_loads.py b/src/fastga/models/load_analysis/wing/aerostructural_loads.py
index 00f83277..b3f28150 100644
--- a/src/fastga/models/load_analysis/wing/aerostructural_loads.py
+++ b/src/fastga/models/load_analysis/wing/aerostructural_loads.py
@@ -173,8 +173,8 @@ def setup(self):
 
     def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
 
-        # STEP 1/XX - DEFINE OR CALCULATE INPUT DATA FOR LOAD COMPUTATION ##############################################
-        ################################################################################################################
+        # STEP 1/XX - DEFINE OR CALCULATE INPUT DATA FOR LOAD COMPUTATION ##########################
+        ############################################################################################
 
         y_vector = inputs["data:aerodynamics:wing:low_speed:Y_vector"]
         y_vector_slip = inputs["data:aerodynamics:slipstream:wing:cruise:prop_on:Y_vector"]
@@ -204,8 +204,9 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         shear_max = 0.0
         rbm_max = 0.0
 
-        # STEP 2/XX - DELETE THE ADDITIONAL ZEROS WE HAD TO PUT TO FIT OPENMDAO AND ADD A POINT AT THE ROOT (Y=0) AND AT
-        # THE VERY TIP (Y=SPAN/2) TO GET THE WHOLE SPAN OF THE WING IN THE INTERPOLATION WE WILL DO LATER ##############
+        # STEP 2/XX - DELETE THE ADDITIONAL ZEROS WE HAD TO PUT TO FIT OPENMDAO AND ADD A POINT
+        # AT THE ROOT (Y=0) AND AT THE VERY TIP (Y=SPAN/2) TO GET THE WHOLE SPAN OF THE WING IN
+        # THE INTERPOLATION WE WILL DO LATER
 
         # We delete the zeros
         y_vector = AerostructuralLoad.delete_additional_zeros(y_vector)
@@ -230,8 +231,8 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         cl_vector_slip = np.append(cl_vector_slip, 0.0)
         chord_vector = np.append(chord_vector, chord_vector[-1])
 
-        # STEP 3/XX - WE COMPUTE THE BASELINE LIFT THAT WE ASSUME WILL SCALE WITH THE LOAD FACTOR, THAT IS WHY WE
-        # COMPUTE It OUT OF THE LOOPS ##################################################################################
+        # STEP 3/XX - WE COMPUTE THE BASELINE LIFT THAT WE ASSUME WILL SCALE WITH THE LOAD
+        # FACTOR, THAT IS WHY WE COMPUTE It OUT OF THE LOOPS
 
         y_vector, _ = self.compute_relief_force(inputs, y_vector_orig, chord_vector, wing_mass, 0.0)
         cl_s = self.compute_cl_s(y_vector_orig, y_vector_orig, y_vector, cl_vector, chord_vector)
@@ -244,16 +245,16 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         weight_shear_diagram = np.full(len(y_vector), 0.0)
         weight_bending_diagram = np.full(len(y_vector), 0.0)
 
-        # STEP 4/XX - WE INITIALIZE THE LOOPS ON THE DIFFERENT SIZING CASE THAT WE DEFINED AND THEN LAUNCH THEM ########
-        ################################################################################################################
+        # STEP 4/XX - WE INITIALIZE THE LOOPS ON THE DIFFERENT SIZING CASE THAT WE DEFINED AND
+        # THEN LAUNCH THEM
 
         mass_array = np.array([mtow, min(mzfw, mtow)])
 
         for mass in mass_array:
 
-            # STEP 4.1/XX - FIRST SUB-STEP IS TO COMPUTE THEN LOAD FACTOR EXPERIENCED BY THE AIRCRAFT AT THE CURRENT
-            # SIZING MASS USING THE FUNCTION WE INHERITED FROM THE ComputeVn CLASS AND THE BASELINE WEIGHT DISTRIBUTION
-            # FOR CURRENT FUEL IN THE WING #############################################################################
+            # STEP 4.1/XX - FIRST SUB-STEP IS TO COMPUTE THEN LOAD FACTOR EXPERIENCED BY THE
+            # AIRCRAFT AT THE CURRENT SIZING MASS USING THE FUNCTION WE INHERITED FROM THE
+            # ComputeVn CLASS AND THE BASELINE WEIGHT DISTRIBUTION FOR CURRENT FUEL IN THE WING
 
             if abs(mtow - mzfw) < 5.0:
                 fuel_mass = 0.0
@@ -279,8 +280,8 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
 
             for load_factor in load_factor_list:
 
-                # STEP 4.2/XX - WE COMPUTE THE REAL CONDITIONS EXPERIENCED IN TERMS OF LIFT AND WEIGHT AND SCALE
-                # THE INITIAL VECTOR ACCORDING TO LOAD FACTOR AND LIFT EQUILIBRIUM
+                # STEP 4.2/XX - WE COMPUTE THE REAL CONDITIONS EXPERIENCED IN TERMS OF LIFT AND
+                # WEIGHT AND SCALE THE INITIAL VECTOR ACCORDING TO LOAD FACTOR AND LIFT EQUILIBRIUM
 
                 cl_wing = 1.05 * (load_factor * mass * 9.81) / (dynamic_pressure * wing_area)
                 cl_s_actual = cl_s * cl_wing / cl_0
@@ -290,9 +291,9 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
                 )
                 weight_array = weight_array_orig * factor_of_safety * load_factor
 
-                # STEP 4.3/XX - WE COMPUTE THE SHEAR AND WEIGHT DIAGRAM WITH THE APPROPRIATE FUNCTION, IDENTIFY THE
-                # MOST EXTREME CONSTRAINTS AND SAVE THE CONDITIONS IN WHICH THEY ARE EXPERIENCED FOR LATER USE IN
-                # THE POST-PROCESSING PHASE ########################################################################
+                # STEP 4.3/XX - WE COMPUTE THE SHEAR AND WEIGHT DIAGRAM WITH THE APPROPRIATE
+                # FUNCTION, IDENTIFY THE MOST EXTREME CONSTRAINTS AND SAVE THE CONDITIONS IN
+                # WHICH THEY ARE EXPERIENCED FOR LATER USE IN THE POST-PROCESSING PHASE
 
                 tot_shear_diagram = AerostructuralLoad.compute_shear_diagram(
                     y_vector, weight_array + lift_section
@@ -323,8 +324,7 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
                     )
                     rbm_max = abs(root_bending_moment)
 
-        # STEP 5/XX - WE ADD ZEROS TO THE RESULTS ARRAYS TO MAKE THEM FIT THE OPENMDAO FORMAT ##########################
-        ################################################################################################################
+        # STEP 5/XX - WE ADD ZEROS TO THE RESULTS ARRAYS TO MAKE THEM FIT THE OPENMDAO FORMAT
 
         additional_zeros = np.zeros(SPAN_MESH_POINT_LOADS - len(y_vector))
         lift_shear_diagram = np.concatenate([lift_shear_diagram, additional_zeros])
@@ -351,9 +351,11 @@ def compute_shear_diagram(y_vector, force_array):
         """
         Function that computes the shear diagram of a given array with linear forces in them
 
-        @param y_vector: an array containing the position of the different station at which the linear forces are given
+        @param y_vector: an array containing the position of the different station at which the
+        linear forces are given
         @param force_array: an array containing the linear forces
-        @return: shear_force_diagram an array representing the shear diagram of the linear forces given in input
+        @return: shear_force_diagram an array representing the shear diagram of the linear forces
+        given in input
         """
 
         # We first create the array to fill
@@ -370,16 +372,19 @@ def compute_bending_moment_diagram(y_vector, force_array):
         """
         Function that computes the root bending diagram of a given array with linear forces in them
 
-        @param y_vector: an array containing the position of the different station at which the linear forces are given
+        @param y_vector: an array containing the position of the different station at which the
+        linear forces are given
         @param force_array: an array containing the linear forces
-        @return: bending_moment_diagram an array representing the root bending diagram of the linear forces given in
+        @return: bending_moment_diagram an array representing the root bending diagram of the
+        linear forces given in
         input
         """
 
         # We first create the array to fill
         bending_moment_diagram = np.zeros(len(y_vector))
 
-        # Each station of the shear diagram is equal to the root bending moment created by all subsequent stations
+        # Each station of the shear diagram is equal to the root bending moment created by all
+        # subsequent stations
         for i in range(len(y_vector)):
             lever_arm = y_vector - y_vector[i]
             bending_moment_diagram[i] = trapz(force_array[i:] * lever_arm[i:], y_vector[i:])
@@ -389,17 +394,21 @@ def compute_bending_moment_diagram(y_vector, force_array):
     @staticmethod
     def compute_cl_s(y_vector_cl_orig, y_vector_chord_orig, y_vector, cl_list, chord_list):
         """
-        Function that computes linear lift on all section of y_vector based on an original cl distribution
-
-        @param y_vector_cl_orig: an array containing the position of the different station at which the original lift
-        distribution was computed, typically a result of OpenVSP or VLM
-        @param y_vector_chord_orig: an array containing the position of the different station at which the chord
-        distribution was computed, typically a result of OpenVSP or VLM
-        @param y_vector: an array containing the position of the different station at which the linear forces are given
+        Function that computes linear lift on all section of y_vector based on an original cl
+        distribution
+
+        @param y_vector_cl_orig: an array containing the position of the different station at which
+        the original lift distribution was computed, typically a result of OpenVSP or VLM
+        @param y_vector_chord_orig: an array containing the position of the different station at
+        which the chord distribution was computed, typically a result of OpenVSP or VLM
+        @param y_vector: an array containing the position of the different station at which the
+        linear forces are given
         @param cl_list: an array containing the original lift coefficient distribution
-        @param chord_list: an array containing the original wing chord length at the different station
-        @return: lift_chord an array representing the linear lift at the different station of y_vector, integrating this
-        vector along the wing span and multiplying it by the dynamic pressure will give you the actual lift distribution
+        @param chord_list: an array containing the original wing chord length at the different
+        station
+        @return: lift_chord an array representing the linear lift at the different station of
+        y_vector, integrating this vector along the wing span and multiplying it by the dynamic
+        pressure will give you the actual lift distribution
         """
 
         # We create the interpolation function
@@ -416,26 +425,29 @@ def compute_cl_s(y_vector_cl_orig, y_vector_chord_orig, y_vector, cl_list, chord
     @staticmethod
     def compute_relief_force(inputs, y_vector, chord_vector, wing_mass, fuel_mass, point_mass=True):
         """
-        Function that computes the baseline weight distribution and modify the y_vector to account for point masses.
-        We chose to represent point masses as linear masses on finite length and to do this we need to modify the
-        y_vector
+        Function that computes the baseline weight distribution and modify the y_vector to
+        account for point masses. We chose to represent point masses as linear masses on finite
+        length and to do this we need to modify the y_vector
 
         @param inputs: inputs parameters defined within FAST-OAD-GA
-        @param y_vector: an array containing the original position of the different station at which the chords are
-        given
+        @param y_vector: an array containing the original position of the different station at which
+        the chords are given
         @param chord_vector: an array containing the chord of the wing at different span station
         @param wing_mass: a float containing the mass of the wing
         @param fuel_mass: a float containing the mass of the fuel
-        @param point_mass: a boolean, if it's FALSE all point mass will be equal to zero used in the post-processing
-        @return: y_vector an array containing the position of the wing span at which the wing mass are sampled
-        @return: weight_array an array containing linear masses of all structural components on the wing
+        @param point_mass: a boolean, if it's FALSE all point mass will be equal to zero used in the
+        post-processing
+        @return: y_vector an array containing the position of the wing span at which the wing mass
+        are sampled
+        @return: weight_array an array containing linear masses of all structural components on the
+        wing
         """
 
-        # STEP 1/XX - DEFINE OR CALCULATE INPUT DATA FOR WEIGHT COMPUTATION ############################################
-        ################################################################################################################
+        # STEP 1/XX - DEFINE OR CALCULATE INPUT DATA FOR WEIGHT COMPUTATION
 
-        # For post-processing we need to be able to compute the linear mass of each component separately, to do this we
-        # chose to render this function able to set each mass to 0 to nullify its influence
+        # For post-processing we need to be able to compute the linear mass of each component
+        # separately, to do this we chose to render this function able to set each mass to 0 to
+        # nullify its influence
         if point_mass:
             tot_engine_mass = inputs["data:weight:propulsion:engine:mass"]
             tot_lg_mass = inputs["data:weight:airframe:landing_gear:main:mass"]
@@ -463,18 +475,17 @@ def compute_relief_force(inputs, y_vector, chord_vector, wing_mass, fuel_mass, p
 
         g = 9.81
 
-        # STEP 2/XX - REARRANGE THE DATA TO FIT ON ONE WING AS WE ASSUME SYMMETRICAL LOADING ###########################
-        ################################################################################################################
+        # STEP 2/XX - REARRANGE THE DATA TO FIT ON ONE WING AS WE ASSUME SYMMETRICAL LOADING
 
         # Computing the mass of the components for one wing
         single_engine_mass = tot_engine_mass / engine_count
         single_lg_mass = tot_lg_mass / 2.0  # We assume 2 MLG
 
-        # STEP 3/XX - AS MENTIONED BEFORE, ADDING POINT MASS WILL ADD Y STATIONS TO THE Y VECTOR SO WE DO THEM
-        # BEFORE ANYTHING ELSE #########################################################################################
+        # STEP 3/XX - AS MENTIONED BEFORE, ADDING POINT MASS WILL ADD Y STATIONS TO THE Y VECTOR
+        # SO WE DO THEM BEFORE ANYTHING ELSE
 
-        # We create the array that will store the "point mass" which we chose to represent as distributed mass over a
-        # small finite interval
+        # We create the array that will store the "point mass" which we chose to represent as
+        # distributed mass over a small finite interval
         point_mass_array = np.zeros(len(y_vector))
 
         # Adding the motor weight
@@ -500,11 +511,12 @@ def compute_relief_force(inputs, y_vector, chord_vector, wing_mass, fuel_mass, p
             y_vector, chord_vector, point_mass_array, y_lg, single_lg_mass, inputs
         )
 
-        # STEP 4/XX - WE CAN NOW ADD THE DISTRIBUTED MASS, I.E THE WING AND THE FUEL. A HARD-CODED VAlUE ENABLE TO
-        # CHANGE FROM ONE MASS DISTRIBUTION TO ANOTHER #################################################################
+        # STEP 4/XX - WE CAN NOW ADD THE DISTRIBUTED MASS, I.E THE WING AND THE FUEL. A
+        # HARD-CODED VAlUE ENABLE TO CHANGE FROM ONE MASS DISTRIBUTION TO ANOTHER
 
-        # We first compute the shape of the distribution regardless of the amplitude which we will later adjust to make
-        # sure that the integration of the mass distribution gives the actual mass
+        # We first compute the shape of the distribution regardless of the amplitude which we
+        # will later adjust to make sure that the integration of the mass distribution gives the
+        # actual mass
         distribution_type = 0.0
 
         if distribution_type == 1.0:
@@ -519,12 +531,12 @@ def compute_relief_force(inputs, y_vector, chord_vector, wing_mass, fuel_mass, p
 
         readjust_fuel = trapz(fuel_weight_distribution, y_vector)
 
-        # We readjust to make sure that the integration of the mass distribution gives the actual mass
+        # We readjust to make sure that the integration of the mass distribution gives the actual
+        # mass
         wing_mass_array = wing_mass * struct_weight_distribution / (2.0 * readjust_struct)
         fuel_mass_array = fuel_mass * fuel_weight_distribution / (2.0 * readjust_fuel)
 
-        # STEP 4/XX - WE CAN NOW ADD ALL THE MASS TOGETHER AND RETURN ALL VALUES #######################################
-        ################################################################################################################
+        # STEP 4/XX - WE CAN NOW ADD ALL THE MASS TOGETHER AND RETURN ALL VALUES
         mass_array = wing_mass_array + fuel_mass_array + point_mass_array
         weight_array = -mass_array * g
 
@@ -537,8 +549,10 @@ def insert_in_sorted_array(array, element):
 
         @param array: a sorted array in which we want to insert an element
         @param element: the element we want to insert in the sorted array
-        @return: final_array a sorted array based on the input array with the argument float inserted in it
-        @return: index the location at which we add to insert the element ot keep the initial array sorted
+        @return: final_array a sorted array based on the input array with the argument float
+        inserted in it
+        @return: index the location at which we add to insert the element ot keep the initial
+        array sorted
         """
 
         tmp_array = np.append(array, element)
@@ -550,12 +564,13 @@ def insert_in_sorted_array(array, element):
     @staticmethod
     def delete_additional_zeros(array, length: int = None):
         """
-        Function that delete the additional zeros we had to add to fit the format imposed by OpenMDAO
+        Function that delete the additional zeros we had to add to fit the format imposed by
+        OpenMDAO
 
         @param array: an array with additional zeros we want to delete
         @param length: if len is specified leave zeros up until the length of the array is len
-        @return: final_array an array containing the same elements of the initial array but with the additional zeros
-        deleted
+        @return: final_array an array containing the same elements of the initial array but with
+        the additional zeros deleted
         """
 
         last_zero = np.amax(np.where(array != 0.0)) + 1
@@ -569,44 +584,45 @@ def delete_additional_zeros(array, length: int = None):
     @staticmethod
     def add_point_mass(y_vector, chord_vector, point_mass_array, y_point_mass, point_mass, inputs):
         """
-        Function that add a point mass to an already created point_mass_array. Modify the y station sampling and chord
-        sampling to account for the additional station added.
+        Function that add a point mass to an already created point_mass_array. Modify the y
+        station sampling and chord sampling to account for the additional station added.
 
-        @param y_vector: the original y_vector which will be modified by adding NB_POINTS_POINT_MASS + 2 points to
-        represent the location of the new point mass
-        @param chord_vector: the original chord vector which will be modified by adding NB_POINTS_POINT_MASS + 2 points
-        to represent the chord at the newly added location
+        @param y_vector: the original y_vector which will be modified by adding
+        NB_POINTS_POINT_MASS + 2 points to represent the location of the new point mass
+        @param chord_vector: the original chord vector which will be modified by adding
+        NB_POINTS_POINT_MASS + 2 points to represent the chord at the newly added location
         @param point_mass_array: the original point mass vector on which we will add the point mass
         @param y_point_mass: the y station of the point mass
         @param point_mass: the value of the mass which we want to add
         @param inputs: inputs parameters defined within FAST-OAD-GA
-        @return: y_vector_new : the new vector contains the y station at which we sample the point mass array with the
-        newly added point mass
+        @return: y_vector_new : the new vector contains the y station at which we sample the point
+        mass array with the newly added point mass
         @return: chord_vector_new : the new vector contains the chord at the new y_station
         @return: point_mass_array_new : the new vector contains the sampled point mass
         """
 
-        # STEP 1/XX - WE EXTRACT THE NECESSARY OUTPUT FROM THE INPUTS ##################################################
-        ################################################################################################################
+        # STEP 1/XX - WE EXTRACT THE NECESSARY OUTPUT FROM THE INPUTS
 
         semi_span = float(inputs["data:geometry:wing:span"]) / 2.0
 
-        # STEP 2/XX - WE CREATE THE INTERPOLATION VECTOR FOR THE CHORD AND THE MASSES. WE ALSO CREATE A VECTOR TO STOCK
-        # WHERE THE MASS ARE ADDED AND HAVE A WAY TO READJUST THE AMPLITUDE ############################################
+        # STEP 2/XX - WE CREATE THE INTERPOLATION VECTOR FOR THE CHORD AND THE MASSES. WE ALSO
+        # CREATE A VECTOR TO STOCK WHERE THE MASS ARE ADDED AND HAVE A WAY TO READJUST THE
+        # AMPLITUDE
 
         fake_point_mass_array = np.zeros(len(point_mass_array))
         present_mass_interp = interp1d(y_vector, point_mass_array)
         present_chord_interp = interp1d(y_vector, chord_vector)
 
-        # STEP 3/XX - WE ALSO STOCK WHERE WE ADD THE Y STATION SINCE IT'LL BE LATER NECESSARY TO READJUST THE AMPLITUDE
-        ################################################################################################################
+        # STEP 3/XX - WE ALSO STOCK WHERE WE ADD THE Y STATION SINCE IT'LL BE LATER NECESSARY TO
+        # READJUST THE AMPLITUDE
 
         interval_len = POINT_MASS_SPAN_RATIO * semi_span / NB_POINTS_POINT_MASS
         nb_point_side = (NB_POINTS_POINT_MASS - 1.0) / 2.0
         y_added = []
 
-        # STEP 4/XX - WE ADD THE NB_POINTS_POINT_MASS AND 2 MORE POINT JUST BEFORE AND AFTER TO GET THE PROPER SQUARE
-        # SHAPE AND NOT A TRAPEZE. WE ALSO ADD THE CORRESPONDING POINT IN THE CHORD VECTOR##############################
+        # STEP 4/XX - WE ADD THE NB_POINTS_POINT_MASS AND 2 MORE POINT JUST BEFORE AND AFTER TO
+        # GET THE PROPER SQUARE SHAPE AND NOT A TRAPEZE. WE ALSO ADD THE CORRESPONDING POINT IN
+        # THE CHORD VECTOR
 
         for i in range(NB_POINTS_POINT_MASS):
             y_current = y_point_mass + (i - nb_point_side) * interval_len
@@ -634,16 +650,17 @@ def add_point_mass(y_vector, chord_vector, point_mass_array, y_point_mass, point
         point_mass_array = np.insert(point_mass_array, index, present_mass_interp(y_max))
         fake_point_mass_array = np.insert(fake_point_mass_array, index, 0.0)
 
-        # STEP 5/XX - WE NOW HAVE THE RIGHT WE JUST NEED TO SCALE IT PROPERLY WHICH IS THE POINT OF THIS STEP ##########
-        ################################################################################################################
+        # STEP 5/XX - WE NOW HAVE THE RIGHT WE JUST NEED TO SCALE IT PROPERLY WHICH IS THE POINT
+        # OF THIS STEP
 
         where_add_mass_grt = np.greater_equal(y_vector, min(y_added))
         where_add_mass_lss = np.less_equal(y_vector, max(y_added))
         where_add_mass = np.logical_and(where_add_mass_grt, where_add_mass_lss)
         where_add_mass_index = np.where(where_add_mass)
 
-        # The fake mass array we use to readjust contains the mass at the station between the first one we added and
-        # the last one which might contains the y station added for other point mass if two are on top of each other
+        # The fake mass array we use to readjust contains the mass at the station between the
+        # first one we added and the last one which might contains the y station added for other
+        # point mass if two are on top of each other
 
         for idx in where_add_mass_index:
             fake_point_mass_array[idx] = 1.0
diff --git a/src/fastga/models/load_analysis/wing/structural_loads.py b/src/fastga/models/load_analysis/wing/structural_loads.py
index 353393d8..221c870d 100644
--- a/src/fastga/models/load_analysis/wing/structural_loads.py
+++ b/src/fastga/models/load_analysis/wing/structural_loads.py
@@ -129,8 +129,7 @@ def setup(self):
 
     def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
 
-        # STEP 1/XX - DEFINE OR CALCULATE INPUT DATA FOR LOAD COMPUTATION ##############################################
-        ################################################################################################################
+        # STEP 1/XX - DEFINE OR CALCULATE INPUT DATA FOR LOAD COMPUTATION
 
         y_vector = inputs["data:aerodynamics:wing:low_speed:Y_vector"]
         chord_vector = inputs["data:aerodynamics:wing:low_speed:chord_vector"]
@@ -144,8 +143,9 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         wing_mass = inputs["data:weight:airframe:wing:mass"]
         fuel_mass = inputs["data:mission:sizing:fuel"]
 
-        # STEP 2/XX - DELETE THE ADDITIONAL ZEROS WE HAD TO PUT TO FIT OPENMDAO AND ADD A POINT AT THE ROOT (Y=0) AND AT
-        # THE VERY TIP (Y=SPAN/2) TO GET THE WHOLE SPAN OF THE WING IN THE INTERPOLATION WE WILL DO LATER ##############
+        # STEP 2/XX - DELETE THE ADDITIONAL ZEROS WE HAD TO PUT TO FIT OPENMDAO AND ADD A POINT
+        # AT THE ROOT (Y=0) AND AT THE VERY TIP (Y=SPAN/2) TO GET THE WHOLE SPAN OF THE WING IN
+        # THE INTERPOLATION WE WILL DO LATER
 
         # Reformat the y_vector array as was done in the aerostructural component
         y_vector = self.delete_additional_zeros(y_vector)
@@ -155,8 +155,8 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         y_vector_orig, _ = self.insert_in_sorted_array(y_vector, semi_span)
         chord_vector_orig = np.append(chord_vector, tip_chord)
 
-        # STEP 3/XX - WE COMPUTE THE BASELINE WEiGHT DISTRIBUTION AND SCALE IT UP ACCORDING TO THE MOST CONSTRAINING
-        # CASE IDENTIFIED IN  THE AEROSTRUCTURAL ANALYSIS ##############################################################
+        # STEP 3/XX - WE COMPUTE THE BASELINE WEiGHT DISTRIBUTION AND SCALE IT UP ACCORDING TO
+        # THE MOST CONSTRAINING CASE IDENTIFIED IN  THE AEROSTRUCTURAL ANALYSIS
 
         y_vector, point_mass_array_orig = self.compute_relief_force(
             inputs, y_vector_orig, chord_vector_orig, 0.0, 0.0
@@ -193,8 +193,8 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
             y_vector, load_factor_shear * fuel_mass_array_orig
         )
 
-        # STEP 4/XX - WE ADD ZEROS AT THE END OF THE RESULT LIFT DISTRIBUTION TO FIT THE FORMAT IMPOSED BY OPENMDAO ####
-        ################################################################################################################
+        # STEP 4/XX - WE ADD ZEROS AT THE END OF THE RESULT LIFT DISTRIBUTION TO FIT THE FORMAT
+        # IMPOSED BY OPENMDAO
 
         point_shear_array = np.concatenate([point_shear_array, additional_zeros])
         wing_shear_array = np.concatenate([wing_shear_array, additional_zeros])
@@ -214,8 +214,8 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
             y_vector, load_factor_rbm * fuel_mass_array_orig
         )
 
-        # STEP 4/XX - WE ADD ZEROS AT THE END OF THE RESULT LIFT DISTRIBUTION TO FIT THE FORMAT IMPOSED BY OPENMDAO ####
-        ################################################################################################################
+        # STEP 4/XX - WE ADD ZEROS AT THE END OF THE RESULT LIFT DISTRIBUTION TO FIT THE FORMAT
+        # IMPOSED BY OPENMDAO
 
         point_root_bending_array = np.concatenate([point_root_bending_array, additional_zeros])
         wing_root_bending_array = np.concatenate([wing_root_bending_array, additional_zeros])
diff --git a/src/fastga/models/loops/update_wing_area_advanced.py b/src/fastga/models/loops/update_wing_area_advanced.py
index 6b71c1dc..5b047714 100644
--- a/src/fastga/models/loops/update_wing_area_advanced.py
+++ b/src/fastga/models/loops/update_wing_area_advanced.py
@@ -20,7 +20,6 @@
 import copy
 
 from scipy.constants import g
-import warnings
 
 from fastoad.module_management.service_registry import RegisterOpenMDAOSystem
 from fastoad.module_management.constants import ModelDomain
@@ -114,14 +113,14 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
 
         if ier != 1:
             _LOGGER.warning(
-                "Could not find a wing area that suits the requirement for fuel inside the wing, using the "
-                "other constraints instead. New value is equal to %f"
-                % max(wing_area_mission, wing_area_approach)
+                "Could not find a wing area that suits the requirement for fuel inside the wing, "
+                "using the other constraints instead. New value is equal to %f",
+                max(wing_area_mission, wing_area_approach),
             )
         else:
             _LOGGER.info(
-                "Looping on wing area with new value equal to %f"
-                % max(wing_area_mission, wing_area_approach)
+                "Looping on wing area with new value equal to %f",
+                max(wing_area_mission, wing_area_approach),
             )
 
         outputs["data:geometry:wing:area"] = max(wing_area_mission, wing_area_approach)
@@ -151,9 +150,10 @@ def compute_wing_area_new(self, wing_area, inputs, fuel_mission):
         y_lg = float(inputs["data:geometry:landing_gear:y"])
         k = float(inputs["settings:geometry:fuel_tanks:depth"])
 
-        # We first have to recompute all the data needed for the tank capacity computation that depends on the wing area
-        # To ensure coherency with the method used, we will use the generate block analysis method on the component that
-        # compute said value and then use a block_analysis on the compute_mfw_advanced component.
+        # We first have to recompute all the data needed for the tank capacity computation that
+        # depends on the wing area To ensure coherency with the method used, we will use the
+        # generate block analysis method on the component that compute said value and then use a
+        # block_analysis on the compute_mfw_advanced component.
 
         # First we need to compute the y positions and the span
 
@@ -240,8 +240,8 @@ def compute_wing_area_new(self, wing_area, inputs, fuel_mission):
         ]
         tip_chord = convert_units(tip_chord_original_value, tip_chord_original_unit, "m")
 
-        # We can now move on to the computation of the mfw for that wing area and then return the difference to solve
-        # for the right wing_area
+        # We can now move on to the computation of the mfw for that wing area and then return the
+        # difference to solve for the right wing_area
 
         var_inputs_compute_mfw = [
             "data:propulsion:IC_engine:fuel_type",
diff --git a/src/fastga/models/loops/update_wing_area_simple.py b/src/fastga/models/loops/update_wing_area_simple.py
index d2164b7b..2257cc79 100644
--- a/src/fastga/models/loops/update_wing_area_simple.py
+++ b/src/fastga/models/loops/update_wing_area_simple.py
@@ -40,7 +40,7 @@ def setup(self):
 
 
 class _UpdateWingAreaSimple(om.ExplicitComponent):
-    """Computation of wing area from needed approach speed and mission fuel"""
+    """Computation of wing area from needed approach speed and mission fuel."""
 
     def setup(self):
         self.add_input("data:mission:sizing:fuel", val=np.nan, units="kg")
@@ -97,8 +97,8 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         wing_area_approach = 2 * mlw * g / (stall_speed ** 2) / (1.225 * max_cl)
 
         _LOGGER.info(
-            "Looping on wing area with new value equal to %f"
-            % max(wing_area_mission, wing_area_approach)
+            "Looping on wing area with new value equal to %f",
+            max(wing_area_mission, wing_area_approach),
         )
 
         outputs["data:geometry:wing:area"] = max(wing_area_mission, wing_area_approach)
diff --git a/src/fastga/models/performances/mission/dynamic_equilibrium.py b/src/fastga/models/performances/mission/dynamic_equilibrium.py
index 1707518d..822f0d03 100644
--- a/src/fastga/models/performances/mission/dynamic_equilibrium.py
+++ b/src/fastga/models/performances/mission/dynamic_equilibrium.py
@@ -49,7 +49,8 @@
 
 class DynamicEquilibrium(om.ExplicitComponent):
     """
-    Compute the derivatives and associated lift-drag-thrust decomposition depending if DP model is included or not
+    Compute the derivatives and associated lift-drag-thrust decomposition depending if DP model
+    is included or not.
     """
 
     def __init__(self, **kwargs):
@@ -134,15 +135,20 @@ def dynamic_equilibrium(
         Method that finds the regulated thrust and aircraft to air angle to obtain dynamic equilibrium
 
         :param inputs: inputs derived from aero and mass models
-        :param gamma: path angle (in rad.) equal to climb rate c=dh/dt over air speed V, sin(gamma)=c/V
+        :param gamma: path angle (in rad.) equal to climb rate c=dh/dt over air speed V,
+        sin(gamma)=c/V
         :param q: dynamic pressure q=1/2*rho*V²
         :param dvx_dt: acceleration linear to air speed
         :param dvz_dt: acceleration perpendicular to air speed
-        :param mass: current mass of the flying aircraft (taking into account propulsion consumption if needed)
-        :param flap_condition: can refer either to "takeoff" or "landing" if high-lift contribution should be considered
+        :param mass: current mass of the flying aircraft (taking into account propulsion consumption
+        if needed)
+        :param flap_condition: can refer either to "takeoff" or "landing" if high-lift contribution
+        should be considered
         :param previous_step: give previous step equilibrium if known to accelerate the calculation
-        :param low_speed: define which aerodynamic models should be used (either low speed or high speed)
-        :param x_cg: x position of the center of gravity of the aircraft, if not given, computed based on fuel in tank
+        :param low_speed: define which aerodynamic models should be used (either low speed or high
+        speed)
+        :param x_cg: x position of the center of gravity of the aircraft, if not given, computed
+        based on fuel in tank
         """
 
         cl_max_clean_htp = inputs["data:aerodynamics:horizontal_tail:low_speed:CL_max_clean"]
@@ -207,8 +213,10 @@ def found_cl_repartition(
         :param mass: current aircraft mass
         :param q: dynamic pressure q=1/2*rho*V²
         :param delta_cm: DP induced cm to be added to the moment equilibrium
-        :param low_speed: define which aerodynamic models should be used (either low speed or high speed)
-        :param x_cg: x_cg position of the aircraft, can be specified. If not specified, computed based on current fuel
+        :param low_speed: define which aerodynamic models should be used (either low speed or high
+        speed)
+        :param x_cg: x_cg position of the aircraft, can be specified. If not specified, computed
+        based on current fuel
         in the aircraft
         """
 
@@ -249,12 +257,18 @@ def found_cl_repartition(
         inv_a = np.linalg.inv(a)
         cl_array = np.dot(inv_a, b)
 
-        # Return equilibrated lift coefficients if low speed maximum clean Cl not exceeded otherwise only cl_wing, 3rd
-        # term is an error flag returned by the function
+        # Return equilibrated lift coefficients if low speed maximum clean Cl not exceeded
+        # otherwise only cl_wing, 3rd term is an error flag returned by the function
         if cl_array[0] < cl_max_clean:
-            return float(cl_array[0]), float(cl_array[1]), False
+            cl_wing_return = float(cl_array[0])
+            cl_htp_return = float(cl_array[1])
+            error = False
         else:
-            return float(mass * g * load_factor / (q * wing_area)), 0.0, True
+            cl_wing_return = float(mass * g * load_factor / (q * wing_area))
+            cl_htp_return = 0.0
+            error = True
+
+        return cl_wing_return, cl_htp_return, error
 
     def save_point(
         self,
@@ -354,9 +368,9 @@ def equation_outer(
         x_cg=None,
     ):
 
-        # Define the system of equations to be solved: load equilibrium along the air x/z axis and
-        # moment equilibrium performed with found_cl_repartition sub-function.
-        # The moment generated by (x_cg_aircraft - x_cg_engine) * T * sin(alpha - alpha_eng) is neglected!
+        # Define the system of equations to be solved: load equilibrium along the air x/z axis
+        # and moment equilibrium performed with found_cl_repartition sub-function. The moment
+        # generated by (x_cg_aircraft - x_cg_engine) * T * sin(alpha - alpha_eng) is neglected!
 
         if low_speed:
             coeff_k_wing = inputs["data:aerodynamics:wing:low_speed:induced_drag_coefficient"]
diff --git a/src/fastga/models/performances/mission/mission.py b/src/fastga/models/performances/mission/mission.py
index 8c928103..89084346 100644
--- a/src/fastga/models/performances/mission/mission.py
+++ b/src/fastga/models/performances/mission/mission.py
@@ -57,6 +57,14 @@ class Mission(om.Group):
 
     """
 
+    def __init__(self, **kwargs):
+        """Defining solvers for mission computation resolution."""
+        super().__init__(**kwargs)
+
+        # Solvers setup
+        self.nonlinear_solver = om.NonlinearBlockGS()
+        self.linear_solver = om.LinearBlockGS()
+
     def initialize(self):
         self.options.declare("propulsion_id", default=None, types=str, allow_none=True)
         self.options.declare("out_file", default="", types=str)
@@ -110,15 +118,11 @@ def setup(self):
         self.add_subsystem("update_fw", UpdateFW(), promotes=["*"])
 
         # Solvers setup
-        self.nonlinear_solver = om.NonlinearBlockGS()
         self.nonlinear_solver.options["debug_print"] = True
-        # self.nonlinear_solver.options["err_on_non_converge"] = True
         self.nonlinear_solver.options["iprint"] = 0
         self.nonlinear_solver.options["maxiter"] = 100
-        # self.nonlinear_solver.options["reraise_child_analysiserror"] = True
         self.nonlinear_solver.options["rtol"] = 1e-2
 
-        self.linear_solver = om.LinearBlockGS()
         # self.linear_solver.options["err_on_non_converge"] = True
         self.linear_solver.options["iprint"] = 0
         self.linear_solver.options["maxiter"] = 10
@@ -198,6 +202,10 @@ def __init__(
     ):
         super().__init__(altitude, delta_t, altitude_in_feet)
         self._calibrated_airspeed = None
+        self._equivalent_airspeed = None
+        self._true_airspeed = None
+        self._mach = None
+        self._unitary_reynolds = None
 
     @property
     def true_airspeed(self) -> Union[float, Sequence[float]]:
@@ -755,7 +763,7 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
 
         # Calculate constant speed (cos(gamma)~1) and corresponding descent angle
         # FIXME: VCAS constant-speed strategy is specific to ICE-propeller configuration, should be an input!
-        atm = _Atmosphere(altitude_t)
+        atm = _Atmosphere(altitude_t, altitude_in_feet=False)
         vs1 = math.sqrt((mass_t * g) / (0.5 * atm.density * wing_area * cl_max_clean))
         v_cas = max(math.sqrt((mass_t * g) / (0.5 * atm.density * wing_area * cl)), 1.3 * vs1)
         atm.calibrated_airspeed = v_cas
diff --git a/src/fastga/models/performances/mission/mission_builder_prep.py b/src/fastga/models/performances/mission/mission_builder_prep.py
index 72441f1f..9660edff 100644
--- a/src/fastga/models/performances/mission/mission_builder_prep.py
+++ b/src/fastga/models/performances/mission/mission_builder_prep.py
@@ -1,4 +1,4 @@
-"""Simple module for prepping the use of the mission builder"""
+"""Simple module for prepping the use of the mission builder."""
 #  This file is part of FAST : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2020  ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -38,8 +38,8 @@ def setup(self):
 
 class _PrepareMissionBuilder(om.ExplicitComponent):
     """
-    Make some simple computation in order to enable the use of the mission builder in FAST-OAD-GA with relevant
-    physic properties
+    Make some simple computation in order to enable the use of the mission builder in FAST-OAD-GA
+    with relevant physic properties
     """
 
     def setup(self):
@@ -67,6 +67,7 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
 
         outputs["data:mission:sizing:cs23:min_climb_speed"] = v_climb_min
 
-        # Based on the ratio between the speed of best L/D and the speed of best endurance according to Gudmundsson
+        # Based on the ratio between the speed of best L/D and the speed of best endurance
+        # according to Gudmundsson
 
         outputs["data:mission:sizing:holding:v_holding"] = 0.75 * inputs["data:TLAR:v_cruise"]
diff --git a/src/fastga/models/performances/mission/takeoff.py b/src/fastga/models/performances/mission/takeoff.py
index cd894a48..d729a000 100644
--- a/src/fastga/models/performances/mission/takeoff.py
+++ b/src/fastga/models/performances/mission/takeoff.py
@@ -27,8 +27,6 @@
 
 # noinspection PyProtectedMember
 from fastoad.module_management._bundle_loader import BundleLoader
-from fastoad.module_management.service_registry import RegisterOpenMDAOSystem
-from fastoad.module_management.constants import ModelDomain
 from fastoad.constants import EngineSetting
 
 from fastga.models.propulsion.fuel_propulsion.base import FuelEngineSet
@@ -211,7 +209,8 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
             else:
                 _LOGGER.critical(
                     "Climb rate is less than %f, adjust weight, propulsion or takeoff "
-                    "configuration" % CLIMB_GRAD_AEO
+                    "configuration",
+                    CLIMB_GRAD_AEO,
                 )
                 raise RuntimeError()
 
diff --git a/src/fastga/models/propulsion/fuel_propulsion/basicIC_engine/basicIC_engine.py b/src/fastga/models/propulsion/fuel_propulsion/basicIC_engine/basicIC_engine.py
index 9b854634..f7aa8bab 100644
--- a/src/fastga/models/propulsion/fuel_propulsion/basicIC_engine/basicIC_engine.py
+++ b/src/fastga/models/propulsion/fuel_propulsion/basicIC_engine/basicIC_engine.py
@@ -188,7 +188,8 @@ def read_map(map_file_path):
         return rpm_vect, pme_vect, pme_limit_vect, sfc_matrix
 
     def compute_flight_points(self, flight_points: FlightPoint):
-        # pylint: disable=too-many-arguments  # they define the trajectory
+        # pylint: disable=too-many-arguments
+        # they define the trajectory
         self.specific_shape = np.shape(flight_points.mach)
         if isinstance(flight_points.mach, float):
             sfc, thrust_rate, thrust = self._compute_flight_points(
diff --git a/src/fastga/models/weight/cg/cg.py b/src/fastga/models/weight/cg/cg.py
index 22233e68..188ece1c 100644
--- a/src/fastga/models/weight/cg/cg.py
+++ b/src/fastga/models/weight/cg/cg.py
@@ -31,6 +31,14 @@
 class CG(om.Group):
     """Model that computes the global center of gravity."""
 
+    def __init__(self, **kwargs):
+        """Defining solvers for cg computation resolution."""
+        super().__init__(**kwargs)
+
+        # Solvers setup
+        self.nonlinear_solver = om.NonlinearBlockGS()
+        self.linear_solver = om.LinearBlockGS()
+
     def initialize(self):
         self.options.declare("propulsion_id", default="", types=str)
 
@@ -49,16 +57,11 @@ def setup(self):
         )
 
         # Solvers setup
-        self.nonlinear_solver = om.NonlinearBlockGS()
         self.nonlinear_solver.options["debug_print"] = True
         self.nonlinear_solver.options["err_on_non_converge"] = True
         self.nonlinear_solver.options["iprint"] = 0
         self.nonlinear_solver.options["maxiter"] = 50
-        # self.nonlinear_solver.options["reraise_child_analysiserror"] = True
-        # self.nonlinear_solver.options["rtol"] = 1e-5
 
-        self.linear_solver = om.LinearBlockGS()
         self.linear_solver.options["err_on_non_converge"] = True
         self.linear_solver.options["iprint"] = 0
         self.linear_solver.options["maxiter"] = 10
-        # self.linear_solver.options["rtol"] = 1e-5
diff --git a/src/fastga/models/weight/cg/cg_components/a_airframe/a1_wing_cg.py b/src/fastga/models/weight/cg/cg_components/a_airframe/a1_wing_cg.py
index e58e1ea3..4ca63a4a 100644
--- a/src/fastga/models/weight/cg/cg_components/a_airframe/a1_wing_cg.py
+++ b/src/fastga/models/weight/cg/cg_components/a_airframe/a1_wing_cg.py
@@ -81,8 +81,9 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
                 chord_reduction = (y_cg - y2_wing) / (y4_wing - y2_wing) * (l2_wing - l4_wing)
 
             chord_at_cg_pos = l2_wing - chord_reduction
-            # In the computation of the mfw, we assume that there is 30% of the chord between front and rear spar,
-            # we will do the same assumption here. It was also assumed that the front spar is at 30 % of the chord
+            # In the computation of the mfw, we assume that there is 30% of the chord between
+            # front and rear spar, we will do the same assumption here. It was also assumed that
+            # the front spar is at 30 % of the chord
             distance_between_spars = 0.30 * chord_at_cg_pos
 
             x_cg_wing_rel = (
diff --git a/src/fastga/models/weight/cg/cg_components/a_airframe/a2_fuselage_cg.py b/src/fastga/models/weight/cg/cg_components/a_airframe/a2_fuselage_cg.py
index 930465be..c37def54 100644
--- a/src/fastga/models/weight/cg/cg_components/a_airframe/a2_fuselage_cg.py
+++ b/src/fastga/models/weight/cg/cg_components/a_airframe/a2_fuselage_cg.py
@@ -60,7 +60,7 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
             x_cg_a2 = lav + 0.35 * (fus_length - lav)
         else:
             _LOGGER.warning(
-                "Propulsion layout %f does not exist, replaced by layout 1!" % prop_layout
+                "Propulsion layout %f does not exist, replaced by layout 1!", prop_layout
             )
             x_cg_a2 = 0.39 * fus_length
 
diff --git a/src/fastga/models/weight/cg/cg_components/b_propulsion/b_cg.py b/src/fastga/models/weight/cg/cg_components/b_propulsion/b_cg.py
index d6907cfd..e7fa53a9 100644
--- a/src/fastga/models/weight/cg/cg_components/b_propulsion/b_cg.py
+++ b/src/fastga/models/weight/cg/cg_components/b_propulsion/b_cg.py
@@ -17,7 +17,6 @@
 
 from .b1_engine_cg import ComputeEngineCG
 from .b2_fuel_lines_cg import ComputeFuelLinesCG
-from .b3_tank_cg import ComputeTankCG
 
 from fastoad.module_management.service_registry import RegisterSubmodel
 
diff --git a/src/fastga/models/weight/cg/cg_components/ratio_aft.py b/src/fastga/models/weight/cg/cg_components/ratio_aft.py
index 24c59181..675c4cdd 100644
--- a/src/fastga/models/weight/cg/cg_components/ratio_aft.py
+++ b/src/fastga/models/weight/cg/cg_components/ratio_aft.py
@@ -39,6 +39,14 @@
     SUBMODEL_AIRCRAFT_X_CG_RATIO, "fastga.submodel.weight.cg.aircraft_empty.x_ratio.legacy"
 )
 class ComputeCGRatioAircraftEmpty(om.Group):
+    def __init__(self, **kwargs):
+        """Defining solvers for aircraft empty cg computation resolution."""
+        super().__init__(**kwargs)
+
+        # Solvers setup
+        self.nonlinear_solver = om.NonlinearBlockGS()
+        self.linear_solver = om.LinearBlockGS()
+
     def setup(self):
         self.add_subsystem(
             "wing_cg", RegisterSubmodel.get_submodel(SUBMODEL_WING_CG), promotes=["*"]
@@ -91,19 +99,14 @@ def setup(self):
         self.add_subsystem("cg_ratio", CGRatio(), promotes=["*"])
 
         # Solvers setup
-        self.nonlinear_solver = om.NonlinearBlockGS()
         self.nonlinear_solver.options["debug_print"] = True
         self.nonlinear_solver.options["err_on_non_converge"] = True
         self.nonlinear_solver.options["iprint"] = 0
         self.nonlinear_solver.options["maxiter"] = 50
-        # self.nonlinear_solver.options["reraise_child_analysiserror"] = True
-        # self.nonlinear_solver.options["rtol"] = 1e-5
 
-        self.linear_solver = om.LinearBlockGS()
         self.linear_solver.options["err_on_non_converge"] = True
         self.linear_solver.options["iprint"] = 0
         self.linear_solver.options["maxiter"] = 10
-        # self.linear_solver.options["rtol"] = 1e-5
 
 
 @RegisterSubmodel(SUBMODEL_AIRCRAFT_X_CG, "fastga.submodel.weight.cg.aircraft_empty.x.legacy")
diff --git a/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_lower_flange.py b/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_lower_flange.py
index 6a207a47..d3166559 100644
--- a/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_lower_flange.py
+++ b/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_lower_flange.py
@@ -109,8 +109,8 @@ def setup(self):
         self.add_input(
             "settings:wing:airfoil:flanges:height_ratio",
             val=0.93,
-            desc="ratio of the height between the two flanges and the the thickness of the airfoil, depends on the "
-            "airfoil",
+            desc="ratio of the height between the two flanges and the the thickness of the "
+            "airfoil, depends on the airfoil",
         )
         self.add_input(
             "settings:materials:aluminium:density",
@@ -161,8 +161,8 @@ def setup(self):
 
     def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         """
-        Component that computes the wing web mass necessary to react to the given linear force vector,
-        according to the methodology developed by Raquel Alonso Castilla.
+        Component that computes the wing web mass necessary to react to the given linear force
+        vector, according to the methodology developed by Raquel Alonso Castilla.
         """
 
         wing_mass = inputs["data:weight:airframe:wing:mass"]
diff --git a/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_ribs_mass.py b/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_ribs_mass.py
index 14545bed..9ea293ea 100644
--- a/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_ribs_mass.py
+++ b/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_ribs_mass.py
@@ -34,7 +34,8 @@ def setup(self):
         self.add_input(
             "settings:wing:airfoil:skin:ka",
             val=0.92,
-            desc="Correction coefficient needed to account for the hypothesis of a rectangular wingbox",
+            desc="Correction coefficient needed to account for the hypothesis of a rectangular "
+            "wingbox",
         )
         self.add_input(
             "settings:wing:airfoil:skin:d_wingbox",
@@ -52,8 +53,8 @@ def setup(self):
 
     def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         """
-        Component that computes the ribs mass necessary to react to the given linear force vector, according to the
-        methodology developed by Raquel Alonso Castilla.
+        Component that computes the ribs mass necessary to react to the given linear force
+        vector, according to the methodology developed by Raquel Alonso Castilla.
         """
         fus_width = inputs["data:geometry:fuselage:maximum_width"]
         fus_height = inputs["data:geometry:fuselage:maximum_height"]
diff --git a/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_secondary_mass.py b/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_secondary_mass.py
index 9eb8a96e..fc175f22 100644
--- a/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_secondary_mass.py
+++ b/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_secondary_mass.py
@@ -26,8 +26,8 @@ def setup(self):
         self.add_input(
             "settings:wing:structure:secondary_mass_ratio",
             val=0.25,
-            desc="Ratio of the mass of the secondary structure and the primary structure (between 0.25 and 0.30 "
-            "according to literature",
+            desc="Ratio of the mass of the secondary structure and the primary structure (between "
+            "0.25 and 0.30 according to literature",
         )
 
         self.add_output("data:weight:airframe:wing:secondary_structure:mass", units="kg")
diff --git a/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_skin_mass.py b/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_skin_mass.py
index cba61ed6..a24048cf 100644
--- a/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_skin_mass.py
+++ b/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_skin_mass.py
@@ -75,7 +75,8 @@ def setup(self):
         self.add_input(
             "settings:wing:airfoil:skin:ka",
             val=0.92,
-            desc="Correction coefficient needed to account for the hypothesis of a rectangular wingbox",
+            desc="Correction coefficient needed to account for the hypothesis of a rectangular "
+            "wingbox",
         )
         self.add_input(
             "settings:wing:airfoil:skin:d_wingbox",
@@ -87,8 +88,8 @@ def setup(self):
 
     def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         """
-        Component that computes the skin mass necessary to react to the given linear force vector, according to the
-        methodology developed by Raquel Alonso Castilla.
+        Component that computes the skin mass necessary to react to the given linear force
+        vector, according to the methodology developed by Raquel Alonso Castilla.
         """
         fus_width = inputs["data:geometry:fuselage:maximum_width"]
         fus_height = inputs["data:geometry:fuselage:maximum_height"]
@@ -107,8 +108,9 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         rho_m = inputs["settings:materials:aluminium:density"]
         shear_modulus = inputs["settings:materials:aluminium:shear_modulus"]
 
-        # There is a slight problem when taper ratio tends to one as it induces a 0 / 0 division that could not be
-        # solved by hand, consequently when taper ratio gets too close to 1. it will be taken as 0.97
+        # There is a slight problem when taper ratio tends to one as it induces a 0 / 0 division
+        # that could not be solved by hand, consequently when taper ratio gets too close to 1. it
+        # will be taken as 0.97
         taper_ratio = inputs["data:geometry:wing:taper_ratio"]
         if taper_ratio > 0.97:
             taper_ratio = 0.97
diff --git a/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_upper_flange.py b/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_upper_flange.py
index 034cd5d9..c3272ecf 100644
--- a/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_upper_flange.py
+++ b/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_upper_flange.py
@@ -109,8 +109,8 @@ def setup(self):
         self.add_input(
             "settings:wing:airfoil:flanges:height_ratio",
             val=0.93,
-            desc="ratio of the height between the two flanges and the the thickness of the airfoil, depends on the "
-            "airfoil",
+            desc="ratio of the height between the two flanges and the the thickness of the "
+            "airfoil, depends on the airfoil",
         )
         self.add_input(
             "settings:materials:aluminium:density",
@@ -161,8 +161,8 @@ def setup(self):
 
     def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         """
-        Component that computes the wing web mass necessary to react to the given linear force vector, according to the
-        methodology developed by Raquel Alonso Castilla.
+        Component that computes the wing web mass necessary to react to the given linear force
+        vector, according to the methodology developed by Raquel Alonso Castilla.
         """
         wing_mass = inputs["data:weight:airframe:wing:mass"]
 
@@ -213,8 +213,9 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
 
         safety_factor = inputs["data:mission:sizing:cs23:safety_factor"]
 
-        # STEP 1 - DELETE THE ADDITIONAL ZEROS WE HAD TO PUT TO FIT OPENMDAO AND ADD A POINT AT THE ROOT (Y=0) AND AT
-        # THE VERY TIP (Y=SPAN/2) TO GET THE WHOLE SPAN OF THE WING IN THE INTERPOLATION WE WILL DO LATER ##############
+        # STEP 1 - DELETE THE ADDITIONAL ZEROS WE HAD TO PUT TO FIT OPENMDAO AND ADD A POINT AT
+        # THE ROOT (Y=0) AND AT THE VERY TIP (Y=SPAN/2) TO GET THE WHOLE SPAN OF THE WING IN THE
+        # INTERPOLATION WE WILL DO LATER
 
         semi_span = wing_span / 2.0
 
diff --git a/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_web_mass.py b/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_web_mass.py
index 83fe0d0e..0a4b5606 100644
--- a/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_web_mass.py
+++ b/src/fastga/models/weight/mass_breakdown/a_airframe/wing_components/compute_web_mass.py
@@ -145,8 +145,8 @@ def setup(self):
 
     def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
         """
-        Component that computes the wing web mass necessary to react to the given linear force vector, according to the
-        methodology developed by Raquel Alonso Castilla.
+        Component that computes the wing web mass necessary to react to the given linear force
+        vector, according to the methodology developed by Raquel Alonso Castilla.
         """
         wing_mass = inputs["data:weight:airframe:wing:mass"]
 
@@ -195,8 +195,9 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
 
         safety_factor = inputs["data:mission:sizing:cs23:safety_factor"]
 
-        # STEP 1 - DELETE THE ADDITIONAL ZEROS WE HAD TO PUT TO FIT OPENMDAO AND ADD A POINT AT THE ROOT (Y=0) AND AT
-        # THE VERY TIP (Y=SPAN/2) TO GET THE WHOLE SPAN OF THE WING IN THE INTERPOLATION WE WILL DO LATER ##############
+        # STEP 1 - DELETE THE ADDITIONAL ZEROS WE HAD TO PUT TO FIT OPENMDAO AND ADD A POINT AT
+        # THE ROOT (Y=0) AND AT THE VERY TIP (Y=SPAN/2) TO GET THE WHOLE SPAN OF THE WING IN THE
+        # INTERPOLATION WE WILL DO LATER
 
         semi_span = wing_span / 2.0
 
diff --git a/src/fastga/models/weight/mass_breakdown/c_systems/sum.py b/src/fastga/models/weight/mass_breakdown/c_systems/sum.py
index 8307d1c7..e602897b 100644
--- a/src/fastga/models/weight/mass_breakdown/c_systems/sum.py
+++ b/src/fastga/models/weight/mass_breakdown/c_systems/sum.py
@@ -27,9 +27,15 @@
 
 @RegisterSubmodel(SUBMODEL_SYSTEMS_MASS, "fastga.submodel.weight.mass.systems.legacy")
 class SystemsWeight(om.Group):
-    """
-    Computes mass of airframe.
-    """
+    """Computes mass of systems."""
+
+    def __init__(self, **kwargs):
+        """Defining solvers for systems weight computation resolution."""
+        super().__init__(**kwargs)
+
+        # Solvers setup
+        self.nonlinear_solver = om.NonlinearBlockGS()
+        self.linear_solver = om.LinearBlockGS()
 
     def setup(self):
         self.add_subsystem(
@@ -70,16 +76,11 @@ def setup(self):
         self.add_subsystem("systems_weight_sum", weight_sum, promotes=["*"])
 
         # Solvers setup
-        self.nonlinear_solver = om.NonlinearBlockGS()
         self.nonlinear_solver.options["debug_print"] = True
         self.nonlinear_solver.options["err_on_non_converge"] = True
         self.nonlinear_solver.options["iprint"] = 0
         self.nonlinear_solver.options["maxiter"] = 50
-        # self.nonlinear_solver.options["reraise_child_analysiserror"] = True
-        # self.nonlinear_solver.options["rtol"] = 1e-3
 
-        self.linear_solver = om.LinearBlockGS()
         self.linear_solver.options["err_on_non_converge"] = True
         self.linear_solver.options["iprint"] = 0
         self.linear_solver.options["maxiter"] = 10
-        # self.linear_solver.options["rtol"] = 1e-3
diff --git a/src/fastga/models/weight/mass_breakdown/mass_breakdown.py b/src/fastga/models/weight/mass_breakdown/mass_breakdown.py
index 45aa9a66..1bc00907 100644
--- a/src/fastga/models/weight/mass_breakdown/mass_breakdown.py
+++ b/src/fastga/models/weight/mass_breakdown/mass_breakdown.py
@@ -49,6 +49,14 @@ class MassBreakdown(om.Group):
                          design payload mass and maximum payload mass must be provided.
     """
 
+    def __init__(self, **kwargs):
+        """Defining solvers for mass breakdown computation resolution."""
+        super().__init__(**kwargs)
+
+        # Solvers setup
+        self.nonlinear_solver = om.NonlinearBlockGS()
+        self.linear_solver = om.LinearBlockGS()
+
     def initialize(self):
         self.options.declare(PAYLOAD_FROM_NPAX, types=bool, default=True)
         self.options.declare("propulsion_id", default="", types=str)
@@ -67,19 +75,14 @@ def setup(self):
         self.add_subsystem("update_mzfw_and_mlw", UpdateMLWandMZFW(), promotes=["*"])
 
         # Solvers setup
-        self.nonlinear_solver = om.NonlinearBlockGS()
         self.nonlinear_solver.options["debug_print"] = True
         self.nonlinear_solver.options["err_on_non_converge"] = True
         self.nonlinear_solver.options["iprint"] = 0
         self.nonlinear_solver.options["maxiter"] = 50
-        # self.nonlinear_solver.options["reraise_child_analysiserror"] = True
-        # self.nonlinear_solver.options["rtol"] = 1e-3
 
-        self.linear_solver = om.LinearBlockGS()
         self.linear_solver.options["err_on_non_converge"] = True
         self.linear_solver.options["iprint"] = 0
         self.linear_solver.options["maxiter"] = 10
-        # self.linear_solver.options["rtol"] = 1e-3
 
 
 @RegisterSubmodel(SUBMODEL_OWE, "fastga.submodel.weight.mass.owe.legacy")
diff --git a/src/fastga/utils/postprocessing/analysis_and_plots.py b/src/fastga/utils/postprocessing/analysis_and_plots.py
index bbb8c5fb..989eb738 100644
--- a/src/fastga/utils/postprocessing/analysis_and_plots.py
+++ b/src/fastga/utils/postprocessing/analysis_and_plots.py
@@ -29,6 +29,8 @@
 from fastoad.io import VariableIO
 from fastoad.openmdao.variables import VariableList
 
+from fastga.models.aerodynamics.components.compute_equilibrated_polar import FIRST_INVALID_COEFF
+
 COLS = plotly.colors.DEFAULT_PLOTLY_COLORS
 
 
@@ -44,8 +46,8 @@ def aircraft_geometry_plot(
     :param name: name to give to the trace added to the figure
     :param fig: existing figure to which add the plot
     :param plot_nacelle: boolean to turn on or off the plotting of the nacelles
-    :param file_formatter: the formatter that defines the format of data file. If not provided, default format will
-                           be assumed.
+    :param file_formatter: the formatter that defines the format of data file. If not provided,
+    default format will be assumed.
     :return: wing plot figure.
     """
     variables = VariableIO(aircraft_file_path, file_formatter).read()
@@ -132,14 +134,18 @@ def aircraft_geometry_plot(
     x_wing = x_wing + wing_25mac_x - 0.25 * wing_mac_length - local_wing_mac_le_x
     x_ht = x_ht + wing_25mac_x + ht_distance_from_wing - local_ht_25mac_x
 
-    # pylint: disable=invalid-name # that's a common naming
+    # pylint: disable=invalid-name
+    # that's a common naming
     x = np.concatenate((x_fuselage, x_wing, x_ht))
-    # pylint: disable=invalid-name # that's a common naming
+    # pylint: disable=invalid-name
+    # that's a common naming
     y = np.concatenate((y_fuselage, y_wing, y_ht))
 
-    # pylint: disable=invalid-name # that's a common naming
+    # pylint: disable=invalid-name
+    # that's a common naming
     y = np.concatenate((-y, y))
-    # pylint: disable=invalid-name # that's a common naming
+    # pylint: disable=invalid-name
+    # that's a common naming
     x = np.concatenate((x, x))
 
     if fig is None:
@@ -268,8 +274,8 @@ def evolution_diagram(
     :param aircraft_file_path: path of data file
     :param name: name to give to the trace added to the figure
     :param fig: existing figure to which add the plot
-    :param file_formatter: the formatter that defines the format of data file. If not provided, default format will
-                           be assumed.
+    :param file_formatter: the formatter that defines the format of data file. If not provided,
+    default format will be assumed.
     :return: V-N plot figure.
     """
     variables = VariableIO(aircraft_file_path, file_formatter).read()
@@ -391,16 +397,16 @@ def cl_wing_diagram(
     file_formatter=None,
 ) -> [go.FigureWidget, go.FigureWidget]:
     """
-    Returns a figure plot of the CL distribution on the semi-wing, and highlights the delta_CL before the added part of
-    the wing or before the reduced part of the wing.
+    Returns a figure plot of the CL distribution on the semi-wing, and highlights the delta_CL
+    before the added part of the wing or before the reduced part of the wing.
 
     :param aircraft_ref_file_path: path of reference aircraft data file
     :param aircraft_mod_file_path: path of modified aircraft data file
     :param prop_on: boolean stating if the rotor is on or off (for single propeller plane)
     :param name_ref: name to give to the trace of the reference aircraft
     :param name_mod: name to give to the trace of the modified aircraft
-    :param file_formatter: the formatter that defines the format of data file. If not provided, default format will
-                           be assumed.
+    :param file_formatter: the formatter that defines the format of data file. If not provided,
+    default format will be assumed.
     :return: Cl distribution figure along the span.
     """
 
@@ -525,8 +531,8 @@ def cg_lateral_diagram(
     :param name: name to give to the trace added to the figure
     :param color: color that we give to the aft, empty and fwd CGs of the aircraft
     :param fig: existing figure to which add the plot
-    :param file_formatter: the formatter that defines the format of data file. If not provided, default format will
-                           be assumed.
+    :param file_formatter: the formatter that defines the format of data file. If not provided,
+    default format will be assumed.
     :return: wing plot figure.
     """
     variables = VariableIO(aircraft_file_path, file_formatter).read()
@@ -989,7 +995,8 @@ def drag_breakdown_diagram(
     nacelle_drag_low_speed = variables["data:aerodynamics:nacelles:low_speed:CD0"].value[0]
     other_drag_low_speed = variables["data:aerodynamics:other:low_speed:CD0"].value[0]
 
-    # CRUD (other undesirable drag). Factor from Gudmundsson book. Introduced in aerodynamics.components.cd0_total.py.
+    # CRUD (other undesirable drag). Factor from Gudmundsson book. Introduced in
+    # aerodynamics.components.cd0_total.py.
     crud_factor = 1.25
 
     if fig is None:
@@ -1104,8 +1111,8 @@ def payload_range(
     :param aircraft_file_path: path of data file
     :param name: name to give to the trace added to the figure
     :param fig: existing figure to which add the plot
-    :param file_formatter: the formatter that defines the format of data file. If not provided, default format will
-                           be assumed.
+    :param file_formatter: the formatter that defines the format of data file. If not provided,
+    default format will be assumed.
     :return: payload range figure.
     """
     variables = VariableIO(aircraft_file_path, file_formatter).read()
@@ -1185,8 +1192,8 @@ def aircraft_polar(
     :param aircraft_file_path: path of data file
     :param name: name to give to the trace added to the figure
     :param fig: existing figure to which add the plot
-    :param file_formatter: the formatter that defines the format of data file. If not provided, default format will
-                           be assumed.
+    :param file_formatter: the formatter that defines the format of data file. If not provided,
+    default format will be assumed.
     :param equilibrated: boolean stating if the polar plotted is the equilibrated one or not
     :return: plane polar figure.
     """
@@ -1194,26 +1201,42 @@ def aircraft_polar(
 
     if equilibrated:
         cl_array_cruise = list(variables["data:aerodynamics:aircraft:cruise:equilibrated:CL"].value)
-        cl_array_cruise = [e for i, e in enumerate(cl_array_cruise) if e != 0]
+        cl_array_cruise = [
+            e for i, e in enumerate(cl_array_cruise) if e != 0 and e < FIRST_INVALID_COEFF
+        ]
         cd_array_cruise = list(variables["data:aerodynamics:aircraft:cruise:equilibrated:CD"].value)
-        cd_array_cruise = [e for i, e in enumerate(cd_array_cruise) if e != 0]
+        cd_array_cruise = [
+            e for i, e in enumerate(cd_array_cruise) if e != 0 and e < FIRST_INVALID_COEFF
+        ]
         cl_array_low_speed = list(
             variables["data:aerodynamics:aircraft:low_speed:equilibrated:CL"].value
         )
-        cl_array_low_speed = [e for i, e in enumerate(cl_array_low_speed) if e != 0]
+        cl_array_low_speed = [
+            e for i, e in enumerate(cl_array_low_speed) if e != 0 and e < FIRST_INVALID_COEFF
+        ]
         cd_array_low_speed = list(
             variables["data:aerodynamics:aircraft:low_speed:equilibrated:CD"].value
         )
-        cd_array_low_speed = [e for i, e in enumerate(cd_array_low_speed) if e != 0]
+        cd_array_low_speed = [
+            e for i, e in enumerate(cd_array_low_speed) if e != 0 and e < FIRST_INVALID_COEFF
+        ]
     else:
         cl_array_cruise = list(variables["data:aerodynamics:aircraft:cruise:CL"].value)
-        cl_array_cruise = [e for i, e in enumerate(cl_array_cruise) if e != 0]
+        cl_array_cruise = [
+            e for i, e in enumerate(cl_array_cruise) if e != 0 and e < FIRST_INVALID_COEFF
+        ]
         cd_array_cruise = list(variables["data:aerodynamics:aircraft:cruise:CD"].value)
-        cd_array_cruise = [e for i, e in enumerate(cd_array_cruise) if e != 0]
+        cd_array_cruise = [
+            e for i, e in enumerate(cd_array_cruise) if e != 0 and e < FIRST_INVALID_COEFF
+        ]
         cl_array_low_speed = list(variables["data:aerodynamics:aircraft:low_speed:CL"].value)
-        cl_array_low_speed = [e for i, e in enumerate(cl_array_low_speed) if e != 0]
+        cl_array_low_speed = [
+            e for i, e in enumerate(cl_array_low_speed) if e != 0 and e < FIRST_INVALID_COEFF
+        ]
         cd_array_low_speed = list(variables["data:aerodynamics:aircraft:low_speed:CD"].value)
-        cd_array_low_speed = [e for i, e in enumerate(cd_array_low_speed) if e != 0]
+        cd_array_low_speed = [
+            e for i, e in enumerate(cd_array_low_speed) if e != 0 and e < FIRST_INVALID_COEFF
+        ]
 
     # Computation of the highest CL/CD ratio which gives the L/D max.
     l_d_max_cruise = max(np.asarray(cl_array_cruise) / np.asarray(cd_array_cruise))


Test patch:
diff --git a/src/fastga/models/aerodynamics/unitary_tests/test_beechcraft_76.py b/src/fastga/models/aerodynamics/unitary_tests/test_beechcraft_76.py
index e5526a98..771d3f17 100644
--- a/src/fastga/models/aerodynamics/unitary_tests/test_beechcraft_76.py
+++ b/src/fastga/models/aerodynamics/unitary_tests/test_beechcraft_76.py
@@ -552,11 +552,11 @@ def test_high_lift():
         XML_FILE,
         delta_cl0_landing=0.5037,
         delta_clmax_landing=0.3613,
-        delta_cm_landing=-0.0680,
+        delta_cm_landing=-0.1552,
         delta_cd_landing=0.005,
         delta_cl0_takeoff=0.1930,
         delta_clmax_takeoff=0.0740,
-        delta_cm_takeoff=-0.0260,
+        delta_cm_takeoff=-0.05949,
         delta_cd_takeoff=0.0004,
         cl_delta_elev=0.5115,
         cd_delta_elev=0.0680,
diff --git a/src/fastga/models/aerodynamics/unitary_tests/test_cirrus_sr22.py b/src/fastga/models/aerodynamics/unitary_tests/test_cirrus_sr22.py
index d995dcbd..1c39c5fa 100644
--- a/src/fastga/models/aerodynamics/unitary_tests/test_cirrus_sr22.py
+++ b/src/fastga/models/aerodynamics/unitary_tests/test_cirrus_sr22.py
@@ -551,11 +551,11 @@ def test_high_lift():
         XML_FILE,
         delta_cl0_landing=0.6167,
         delta_clmax_landing=0.4194,
-        delta_cm_landing=-0.1049,
+        delta_cm_landing=-0.1729,
         delta_cd_landing=0.0136,
         delta_cl0_takeoff=0.2300,
         delta_clmax_takeoff=0.0858,
-        delta_cm_takeoff=-0.0391,
+        delta_cm_takeoff=-0.0645,
         delta_cd_takeoff=0.0011,
         cl_delta_elev=0.4686,
         cd_delta_elev=0.06226,
diff --git a/src/fastga/models/performances/unitary_tests/test_beechcraft_76.py b/src/fastga/models/performances/unitary_tests/test_beechcraft_76.py
index bc1e6324..0796187f 100644
--- a/src/fastga/models/performances/unitary_tests/test_beechcraft_76.py
+++ b/src/fastga/models/performances/unitary_tests/test_beechcraft_76.py
@@ -250,11 +250,11 @@ def test_compute_descent():
     # Run problem and check obtained value(s) is/(are) correct
     problem = run_system(group, ivc)
     fuel_mass = problem.get_val("data:mission:sizing:main_route:descent:fuel", units="kg")
-    assert fuel_mass == pytest.approx(0.81, abs=1e-2)
+    assert fuel_mass == pytest.approx(1.12, abs=1e-2)
     distance = (
         problem.get_val("data:mission:sizing:main_route:descent:distance", units="m") / 1000
     )  # conversion to km
-    assert distance == pytest.approx(84, abs=1)
+    assert distance == pytest.approx(91, abs=1)
     duration = problem.get_val("data:mission:sizing:main_route:descent:duration", units="min")
     assert duration == pytest.approx(27, abs=1)
 
diff --git a/src/fastga/models/performances/unitary_tests/test_cirrus_sr22.py b/src/fastga/models/performances/unitary_tests/test_cirrus_sr22.py
index fba98468..6e857ffc 100644
--- a/src/fastga/models/performances/unitary_tests/test_cirrus_sr22.py
+++ b/src/fastga/models/performances/unitary_tests/test_cirrus_sr22.py
@@ -255,11 +255,11 @@ def test_compute_descent():
     # Run problem and check obtained value(s) is/(are) correct
     problem = run_system(group, ivc)
     fuel_mass = problem.get_val("data:mission:sizing:main_route:descent:fuel", units="kg")
-    assert fuel_mass == pytest.approx(0.61, abs=1e-2)
+    assert fuel_mass == pytest.approx(0.83, abs=1e-2)
     distance = (
         problem.get_val("data:mission:sizing:main_route:descent:distance", units="m") / 1000
     )  # conversion to km
-    assert distance == pytest.approx(80, abs=1)
+    assert distance == pytest.approx(87, abs=1)
     duration = problem.get_val("data:mission:sizing:main_route:descent:duration", units="min")
     assert duration == pytest.approx(27, abs=1)
 
@@ -300,7 +300,7 @@ def test_payload_range():
     payload_result = np.array([420, 420, 355, 310.697, 0])
     assert np.max(np.abs(payload_array - payload_result)) <= 1e-1
     range_array = problem.get_val("data:payload_range:range_array", units="NM")
-    range_result = np.array([0.0, 1169.73, 1632.40, 1955.12, 2210.96])
+    range_result = np.array([0.0, 1171.88, 1634.65, 1957.44, 2213.56])
     assert np.max(np.abs(range_array - range_result)) <= 1e-1
     specific_range_array = problem.get_val("data:payload_range:specific_range_array", units="NM/kg")
     specific_range_result = np.array([0.0, 6.24, 6.47, 6.58, 7.45])