diff --git a/src/thermohl/power/__init__.py b/src/thermohl/power/__init__.py
index 3428bf56..8391ba2c 100644
--- a/src/thermohl/power/__init__.py
+++ b/src/thermohl/power/__init__.py
@@ -8,7 +8,7 @@
 from .radiative_cooling import RadiativeCoolingBase
 
 from .power_term import PowerTerm
-from .solar_heating import _SRad, SolarHeatingBase
+from .solar_heating import _SRad, SolarHeatingBase, FixedSolarIrradianceSolarHeating
 
 
 __all__ = [
@@ -16,4 +16,5 @@
     "PowerTerm",
     "_SRad",
     "SolarHeatingBase",
+    "FixedSolarIrradianceSolarHeating",
 ]
diff --git a/src/thermohl/power/cigre/convective_cooling.py b/src/thermohl/power/cigre/convective_cooling.py
index ae95d0ed..563db3c0 100644
--- a/src/thermohl/power/cigre/convective_cooling.py
+++ b/src/thermohl/power/cigre/convective_cooling.py
@@ -11,15 +11,16 @@
 import numpy as np
 
 from thermohl import floatArrayLike
-from thermohl.power import PowerTerm
-from thermohl.power.convective_cooling import compute_wind_attack_angle
+from thermohl.power.convective_cooling import (
+    ConvectiveCoolingBase,
+)
 from thermohl.power.cigre import Air
 
 
 logger = logging.getLogger(__name__)
 
 
-class ConvectiveCooling(PowerTerm):
+class ConvectiveCooling(ConvectiveCoolingBase):
     """Convective cooling term."""
 
     def __init__(
@@ -30,8 +31,8 @@ def __init__(
         wind_speed: floatArrayLike,
         outer_diameter: floatArrayLike,
         roughness_ratio: floatArrayLike,
-        wind_azimuth: floatArrayLike = None,
-        wind_attack_angle: floatArrayLike = None,
+        wind_azimuth: floatArrayLike | None = None,
+        wind_attack_angle: floatArrayLike | None = None,
         g: float = 9.81,
         **kwargs: Any,
     ):
@@ -44,12 +45,16 @@ def __init__(
             cable_azimuth (float | numpy.ndarray): Azimuth (deg).
             ambient_temperature (float | numpy.ndarray): Ambient temperature (°C).
             wind_speed (float | numpy.ndarray): Wind speed (m·s⁻¹).
-            wind_azimuth (float | numpy.ndarray): wind azimuth regarding north (deg).
+            wind_azimuth (float | numpy.ndarray | None): wind azimuth regarding north (deg).
+            wind_attack_angle (float | numpy.ndarray | None): wind attack angle (rad).
             outer_diameter (float | numpy.ndarray): External diameter (m).
             roughness_ratio (float | numpy.ndarray): Cable roughness (—).
             g (float, optional): Gravitational acceleration (m·s⁻²). The default is 9.81.
 
         """
+        self._check_arguments(wind_azimuth, wind_attack_angle)
+        self._set_wind_attack_angle(cable_azimuth, wind_azimuth, wind_attack_angle)
+
         self.altitude = altitude
         self.ambient_temp = ambient_temperature
         self.wind_speed = wind_speed
@@ -57,20 +62,7 @@ def __init__(
         self.roughness_ratio = roughness_ratio
         self.gravity = g
 
-        if wind_attack_angle is None and wind_azimuth is None:
-            raise ValueError("Must provide either wind_attack_angle or wind_azimuth.")
-        if wind_attack_angle is not None and wind_azimuth is not None:
-            logger.warning(
-                "both wind_attack_angle and wind_azimuth are provided. wind_azimuth will be ignored."
-            )
-        if wind_attack_angle is not None:
-            self.wind_attack_angle = wind_attack_angle
-        else:
-            self.wind_attack_angle = compute_wind_attack_angle(
-                cable_azimuth, wind_azimuth
-            )
-
-    def _nu_forced(
+    def _nusselt_forced(
         self, film_temperature: floatArrayLike, kinematic_viscosity: floatArrayLike
     ) -> floatArrayLike:
         """
@@ -133,7 +125,7 @@ def _nu_forced(
             B1 * reynolds**n
         )
 
-    def _nu_natural(
+    def _nusselt_natural(
         self,
         film_temperature: floatArrayLike,
         temperature_delta: floatArrayLike,
@@ -192,8 +184,8 @@ def value(self, conductor_temperature: floatArrayLike) -> floatArrayLike:
         # nu[nu < 1.0E-06] = 1.0E-06
         thermal_conductivity = Air.thermal_conductivity(film_temperature)
         # lm[lm < 0.01] = 0.01
-        nusselt_forced = self._nu_forced(film_temperature, kinematic_viscosity)
-        nusselt_natural = self._nu_natural(
+        nusselt_forced = self._nusselt_forced(film_temperature, kinematic_viscosity)
+        nusselt_natural = self._nusselt_natural(
             film_temperature, temperature_delta, kinematic_viscosity
         )
         return (
diff --git a/src/thermohl/power/convective_cooling.py b/src/thermohl/power/convective_cooling.py
index 65c3159a..784296f9 100644
--- a/src/thermohl/power/convective_cooling.py
+++ b/src/thermohl/power/convective_cooling.py
@@ -46,32 +46,77 @@ def __init__(
         air_density: Callable[[floatArrayLike, floatArrayLike], floatArrayLike],
         dynamic_viscosity: Callable[[floatArrayLike], floatArrayLike],
         thermal_conductivity: Callable[[floatArrayLike], floatArrayLike],
-        wind_azimuth: floatArrayLike = None,
-        wind_attack_angle: floatArrayLike = None,
+        wind_azimuth: floatArrayLike | None = None,
+        wind_attack_angle: floatArrayLike | None = None,
         **kwargs: Any,
     ):
+        self._check_arguments(wind_azimuth, wind_attack_angle)
+        self._set_wind_attack_angle(cable_azimuth, wind_azimuth, wind_attack_angle)
+
         self.altitude = altitude
         self.ambient_temp = ambient_temperature
         self.wind_speed = wind_speed
 
-        if wind_attack_angle is None and wind_azimuth is None:
+        self.outer_diameter = outer_diameter
+
+        self.air_density = air_density
+        self.dynamic_viscosity = dynamic_viscosity
+        self.thermal_conductivity = thermal_conductivity
+
+    @classmethod
+    def _check_arguments(
+        cls,
+        wind_azimuth: floatArrayLike | None,
+        wind_attack_angle: floatArrayLike | None,
+    ) -> None:
+        if (wind_attack_angle is None or np.isnan(wind_attack_angle).any()) and (
+            wind_azimuth is None or np.isnan(wind_azimuth).any()
+        ):
             raise ValueError("Must provide either wind_attack_angle or wind_azimuth.")
-        if wind_attack_angle is not None and wind_azimuth is not None:
+
+        if (
+            wind_attack_angle is not None
+            and not np.isnan(wind_attack_angle).all()
+            and wind_azimuth is not None
+        ):
             logger.warning(
-                "both wind_attack_angle and wind_azimuth are provided. wind_azimuth will be ignored."
+                "Both wind_attack_angle and wind_azimuth are provided. wind_azimuth will be ignored."
             )
-        if wind_attack_angle is not None:
+
+    def _set_wind_attack_angle(
+        self,
+        cable_azimuth: floatArrayLike,
+        wind_azimuth: floatArrayLike,
+        wind_attack_angle: floatArrayLike,
+    ) -> None:
+        # Compute missing wind attack angles
+        if isinstance(wind_attack_angle, np.ndarray) and wind_attack_angle.ndim > 0:
+            if (
+                not isinstance(cable_azimuth, np.ndarray)
+                or cable_azimuth.shape != wind_attack_angle.shape
+            ):
+                raise ValueError(
+                    "If wind_attack_angle is an array, cable_azimuth must be an array of the same shape."
+                )
+            if (
+                not isinstance(wind_azimuth, np.ndarray)
+                or wind_azimuth.shape != wind_attack_angle.shape
+            ):
+                raise ValueError(
+                    "If wind_attack_angle is an array, wind_azimuth must be an array of the same shape."
+                )
+            mask = np.isnan(wind_attack_angle)
+            if np.any(mask):
+                wind_attack_angle[mask] = compute_wind_attack_angle(
+                    cable_azimuth[mask], wind_azimuth[mask]
+                )
             self.wind_attack_angle = wind_attack_angle
-        else:
+        elif wind_attack_angle is None or np.isnan(wind_attack_angle):
             self.wind_attack_angle = compute_wind_attack_angle(
                 cable_azimuth, wind_azimuth
             )
-
-        self.outer_diameter = outer_diameter
-
-        self.air_density = air_density
-        self.dynamic_viscosity = dynamic_viscosity
-        self.thermal_conductivity = thermal_conductivity
+        else:
+            self.wind_attack_angle = wind_attack_angle
 
     def _value_forced(
         self,
diff --git a/src/thermohl/power/solar_heating.py b/src/thermohl/power/solar_heating.py
index af0d953f..a6f3f9d4 100644
--- a/src/thermohl/power/solar_heating.py
+++ b/src/thermohl/power/solar_heating.py
@@ -147,3 +147,21 @@ def derivative(self, conductor_temperature: floatArrayLike) -> floatArrayLike:
         :return: Derivative of solar heating.
         """
         return np.zeros_like(conductor_temperature)
+
+
+class FixedSolarIrradianceSolarHeating(SolarHeatingBase):
+    """Solar heating term with fixed solar irradiance.
+
+    This class only computes the solar heating (power) and its
+    derivative based on the provided solar irradiance.
+    """
+
+    def __init__(
+        self,
+        outer_diameter: floatArrayLike,
+        solar_absorptivity: floatArrayLike,
+        solar_irradiance: floatArrayLike,
+    ):
+        self.outer_diameter = outer_diameter
+        self.solar_absorptivity = solar_absorptivity
+        self.solar_irradiance = solar_irradiance
diff --git a/src/thermohl/solver/parameters.py b/src/thermohl/solver/parameters.py
index be7f970a..03a86295 100644
--- a/src/thermohl/solver/parameters.py
+++ b/src/thermohl/solver/parameters.py
@@ -51,6 +51,7 @@ def _set_default_values(self) -> None:
         self.precipitation_rate = 0.0  # rain precipitation rate (m.s**-1)
         self.wind_speed = 0.0  # wind speed (m.s**-1)
         self.wind_azimuth = 90.0  # wind_azimuth (deg, regarding north)
+        self.wind_attack_angle = np.nan  # wind attack angle (rad)
         self.nebulosity = np.nan  # nebulosity (1)
         self.albedo = 0.15  # albedo (1)
         # coefficient for air pollution from 0 (clean) to 1 (polluted)
diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py
index 769ca840..ff19718d 100644
--- a/src/thermohl/solver/slv1t.py
+++ b/src/thermohl/solver/slv1t.py
@@ -5,16 +5,22 @@
 # file, You can obtain one at http://mozilla.org/MPL/2.0/.
 # SPDX-License-Identifier: MPL-2.0
 
+import logging
 import numbers
 from typing import Optional
 
 import numpy as np
 
 from thermohl import floatArrayLike, floatArray
+from thermohl.power import FixedSolarIrradianceSolarHeating
 from thermohl.solver.solver import Solver as Solver_, get_time_changing_parameters
 from thermohl.solver.parameters import DEFAULT_PARAMETERS as default
 from thermohl.solver.entities import PowerType, VariableType
 from thermohl.utils import bisect_v
+from thermohl.utils import quasi_newton
+
+
+logger = logging.getLogger(__name__)
 
 
 class Solver1T(Solver_):
@@ -220,6 +226,112 @@ def fun(i: floatArray) -> floatArrayLike:
             return_power,
         )
 
-        result = self._add_input_data_to_result(result)
+        return self._add_input_data_to_result(result)
 
-        return result
+    def _set_default_reduced_intensity_args(
+        self,
+        ambient_temperature: Optional[floatArrayLike],
+        wind_speed: Optional[floatArrayLike],
+        solar_irradiance: Optional[floatArrayLike],
+    ):
+        if ambient_temperature is None:
+            logger.warning(
+                "ambient_temperature is not set. Using default value of 30 °C."
+            )
+            ambient_temperature = 30.0
+        self.args.ambient_temperature = ambient_temperature
+
+        if wind_speed is None:
+            logger.warning("wind_speed is not set. Using default value of 0.6 m/s.")
+            wind_speed = 0.6
+        self.args.wind_speed = wind_speed
+
+        if solar_irradiance is None:
+            logger.warning(
+                "solar_irradiance is not set. Using default value of 600 W/m²."
+            )
+            solar_irradiance = 600.0
+        return solar_irradiance
+
+    def reduced_intensity(
+        self,
+        measured_temperature_difference: floatArrayLike,
+        measured_intensity: floatArrayLike,
+        ambient_temperature: Optional[floatArrayLike] = None,
+        wind_speed: Optional[floatArrayLike] = None,
+        solar_irradiance: Optional[floatArrayLike] = None,
+        max_conductor_temperature: Optional[floatArrayLike] = None,
+    ) -> floatArrayLike:
+        """
+        Compute the reduced intensity limit for a given measured temperature difference
+        betwwen the sound cable and a hotspot on the junction between a cable
+        and a faulty sleeve.
+
+        Args:
+            measured_temperature_difference (float | np.ndarray): The measured temperature difference between the cable surface and the sleeve.
+            measured_intensity (float | np.ndarray): The measuredintensity at which the temperature difference was measured.
+            ambient_temperature (Optional[float | np.ndarray]): The ambient temperature. Default is 30.
+            wind_speed (Optional[float | np.ndarray]): The wind speed. Default is 0.6.
+            solar_irradiance (Optional[float | np.ndarray]): The measured solar irradiance. Default is 600.
+            max_conductor_temperature (Optional[float | np.ndarray]): The maximum conductor temperature. Default is 100.
+        """
+        # Save elements that will be modified to do the computations
+        # so as to be able to restore afterwards
+        saved_solver_transit = self.args.transit
+        saved_solver_ambient_temperature = self.args.ambient_temperature
+        saved_solver_wind_speed = self.args.wind_speed
+        saved_solver_wind_attack_angle = self.args.wind_attack_angle
+
+        saved_solar_heating = self.solar_heating
+
+        # Set args default values for reduced intensity computation.
+        # These differ from those used for the other computations.
+        solar_irradiance = self._set_default_reduced_intensity_args(
+            ambient_temperature, wind_speed, solar_irradiance
+        )
+
+        # Set default value for max_conductor_temperature if not provided.
+        if max_conductor_temperature is None:
+            max_conductor_temperature = np.full_like(measured_intensity, 100.0)
+
+        self.args.wind_attack_angle = 90.0
+        self.convective_cooling.__init__(**self.args.__dict__)
+
+        self.solar_heating = FixedSolarIrradianceSolarHeating(
+            outer_diameter=self.args.outer_diameter,
+            solar_absorptivity=self.args.solar_absorptivity,
+            solar_irradiance=solar_irradiance,
+        )
+
+        def conductor_temperature(transit):
+            self.args.transit = transit
+            self.joule_heating.__init__(**self.args.__dict__)
+            return self.steady_temperature()[VariableType.TEMPERATURE][0]
+
+        def temperature_difference(transit):
+            return measured_temperature_difference * (
+                (transit / measured_intensity) ** 2
+            )
+
+        def sleeve_temperature(transit):
+            return conductor_temperature(transit) + temperature_difference(transit)
+
+        def f(transit):
+            return sleeve_temperature(transit) - max_conductor_temperature
+
+        x0 = np.full_like(measured_intensity, 100.0)
+
+        reduced_intensity = quasi_newton(f, x0=x0)
+
+        # Restore previous elements
+        self.args.transit = saved_solver_transit
+        # Update joule heating with restored transit
+        self.joule_heating.__init__(**self.args.__dict__)
+        self.args.ambient_temperature = saved_solver_ambient_temperature
+        self.args.wind_speed = saved_solver_wind_speed
+        self.args.wind_attack_angle = saved_solver_wind_attack_angle
+        # Update convective cooling with restored wind_attack_angle
+        self.convective_cooling.__init__(**self.args.__dict__)
+        self.solar_heating = saved_solar_heating
+
+        return reduced_intensity
diff --git a/src/thermohl/utils.py b/src/thermohl/utils.py
index 84c66d08..86cb2d33 100644
--- a/src/thermohl/utils.py
+++ b/src/thermohl/utils.py
@@ -13,10 +13,14 @@
 import os
 from functools import wraps
 from importlib.util import find_spec
+import warnings
+from typing import Callable
 
 import numpy as np
 import yaml
 
+from thermohl import floatArrayLike
+
 
 logger = logging.getLogger(__name__)
 
@@ -188,6 +192,132 @@ def bisect_v(
     return midpoint, abs_error
 
 
+def _array_quasi_newton(
+    func: Callable[[np.ndarray], np.ndarray], x0: np.ndarray, tol: float, maxiter: int
+) -> floatArrayLike:
+    """
+    A vectorized version of secant method for arrays.
+
+    Do not use this method directly. This method is called from `quasi_newton`
+    when ``np.size(x0) > 1`` is ``True``.
+
+    Heavily inspired by the implementation of the SciPy library.
+    """
+    # Explicitly copy `x0` as `p` will be modified inplace, but the
+    # user's array should not be altered.
+    p = np.array(x0, copy=True)
+
+    failures = np.ones_like(p, dtype=bool)
+    nz_der = np.ones_like(failures)
+
+    dx = np.finfo(float).eps ** 0.33
+    p1 = p * (1 + dx) + np.where(p >= 0, dx, -dx)
+    q0 = np.asarray(func(p))
+    q1 = np.asarray(func(p1))
+    active = np.ones_like(p, dtype=bool)
+    for _ in range(maxiter):
+        nz_der = q1 != q0
+        # stop iterating if all derivatives are zero
+        if not nz_der.any():
+            p = (p1 + p) / 2.0
+            break
+        # Secant Step
+        dp = (q1 * (p1 - p))[nz_der] / (q1 - q0)[nz_der]
+        # only update nonzero derivatives
+        p = np.asarray(p, dtype=np.result_type(p, p1, dp, np.float64))
+        p[nz_der] = p1[nz_der] - dp
+        active_zero_der = ~nz_der & active
+        p[active_zero_der] = (p1 + p)[active_zero_der] / 2.0
+        active &= nz_der  # don't assign zero derivatives again
+        failures[nz_der] = np.abs(dp) >= tol  # not yet converged
+        # stop iterating if there aren't any failures, not incl zero der
+        if not failures[nz_der].any():
+            break
+        p1, p = p, p1
+        q0 = q1
+        q1 = np.asarray(func(p1))
+
+    zero_der = ~nz_der & failures  # don't include converged with zero-ders
+    if zero_der.any():
+        nonzero_dp = p1 != p
+        # non-zero dp, but infinite newton step
+        zero_der_nz_dp = zero_der & nonzero_dp
+        if zero_der_nz_dp.any():
+            rms = np.sqrt(sum((p1[zero_der_nz_dp] - p[zero_der_nz_dp]) ** 2))
+            warnings.warn(f"RMS of {rms:g} reached", RuntimeWarning, stacklevel=3)
+    elif failures.any():
+        all_or_some = "all" if failures.all() else "some"
+        msg = f"{all_or_some:s} failed to converge after {maxiter:d} iterations"
+        if failures.all():
+            raise RuntimeError(msg)
+        warnings.warn(msg, RuntimeWarning, stacklevel=3)
+
+    return p
+
+
+def _check_quasi_newton_arguments(tol: float, maxiter: int, rtol: float) -> None:
+    if tol <= 0:
+        raise ValueError(f"tol too small ({tol:g} <= 0)")
+    if maxiter < 1:
+        raise ValueError("maxiter must be greater than 0")
+    if rtol < 0:
+        raise ValueError(f"rtol too small ({rtol:g} < 0)")
+
+
+def quasi_newton(  # NOSONAR(S3776)
+    func: Callable[[floatArrayLike], floatArrayLike],
+    x0: floatArrayLike,
+    tol: float = 1.48e-8,
+    maxiter: int = 50,
+    rtol: float = 0.0,
+) -> floatArrayLike:
+    """Find the zero of a function using the quasi-Newton (secant) method.
+
+    Heavily inspired by the implementation of optimize.newton in the SciPy library.
+    """
+    _check_quasi_newton_arguments(tol, maxiter, rtol)
+    if np.size(x0) > 1:
+        return _array_quasi_newton(func, x0, tol, maxiter)
+
+    # Convert to float (don't use float(x0); this works also for complex x0)
+    # Use np.asarray because we want x0 to be a numpy object, not a Python
+    # object. e.g. np.complex(1+1j) > 0 is possible, but (1 + 1j) > 0 raises
+    # a TypeError
+    x0 = np.asarray(x0)[()] * 1.0
+    p0 = x0
+
+    eps = 1e-4
+    p1 = x0 * (1 + eps)
+    p1 += eps if p1 >= 0 else -eps
+    q0 = func(p0)
+    q1 = func(p1)
+    if abs(q1) < abs(q0):
+        p0, p1, q0, q1 = p1, p0, q1, q0
+    for itr in range(maxiter):
+        if q1 == q0:
+            if p1 != p0:
+                msg = f"Tolerance of {p1 - p0} reached."
+                msg += (
+                    f" Failed to converge after {itr + 1} iterations,"
+                    f" value is {p1}."
+                )
+                raise RuntimeError(msg)
+            return (p1 + p0) / 2.0
+        else:
+            if abs(q1) > abs(q0):
+                p = (-q0 / q1 * p1 + p0) / (1 - q0 / q1)
+            else:
+                p = (-q1 / q0 * p0 + p1) / (1 - q1 / q0)
+        if np.isclose(p, p1, rtol=rtol, atol=tol):
+            return p
+        p0, q0 = p1, q1
+        p1 = p
+        q1 = func(p1)
+
+    msg = f"Failed to converge after {itr + 1} iterations, value is {p}."
+    raise RuntimeError(msg)
+
+
 # In agreement with Eurobios, this function has been retrieved from the pyntb library,
 # in order to remove the external dependency on this library.
 # In this library, this function was initially developed under the name qnewt2d_v
diff --git a/test/unit/power/cigre/test_power_cigre_convective_cooling.py b/test/unit/power/cigre/test_power_cigre_convective_cooling.py
index 3a746f55..c3401112 100644
--- a/test/unit/power/cigre/test_power_cigre_convective_cooling.py
+++ b/test/unit/power/cigre/test_power_cigre_convective_cooling.py
@@ -45,12 +45,12 @@
     conv_cool_instances,
     ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"],
 )
-def test_nu_forced_float_value(convective_cooling, expected_type):
+def test_nusselt_forced_float_value(convective_cooling, expected_type):
     Tf = 30.0
     nu = 1.5e-5
     expected_result = 17.3425
 
-    result = convective_cooling._nu_forced(Tf, nu)
+    result = convective_cooling._nusselt_forced(Tf, nu)
 
     np.testing.assert_allclose(result, expected_result, rtol=1e-5)
 
@@ -60,12 +60,12 @@ def test_nu_forced_float_value(convective_cooling, expected_type):
     conv_cool_instances,
     ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"],
 )
-def test_nu_forced_array_single_value(convective_cooling, expected_type):
+def test_nusselt_forced_array_single_value(convective_cooling, expected_type):
     Tf = np.array([30.0])
     nu = np.array([1.5e-5])
     expected_result = np.array([17.3425])
 
-    result = convective_cooling._nu_forced(Tf, nu)
+    result = convective_cooling._nusselt_forced(Tf, nu)
 
     np.testing.assert_allclose(result, expected_result, rtol=1e-5)
 
@@ -75,12 +75,12 @@ def test_nu_forced_array_single_value(convective_cooling, expected_type):
     conv_cool_instances,
     ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"],
 )
-def test_nu_forced_array_values(convective_cooling, expected_type):
+def test_nusselt_forced_array_values(convective_cooling, expected_type):
     Tf = np.array([30.0, 35.0])
     nu = np.array([1.5e-5, 1.6e-5])
     expected_result = np.array([17.3425, 16.6482])
 
-    result = convective_cooling._nu_forced(Tf, nu)
+    result = convective_cooling._nusselt_forced(Tf, nu)
 
     np.testing.assert_allclose(result, expected_result, rtol=1e-5)
 
@@ -90,14 +90,14 @@ def test_nu_forced_array_values(convective_cooling, expected_type):
     conv_cool_instances,
     ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"],
 )
-def test_nu_forced_boundary_conditions(convective_cooling, expected_type):
+def test_nusselt_forced_boundary_conditions(convective_cooling, expected_type):
     Tf = np.array([30.0])
     nu = np.array([1.5e-5])
     convective_cooling.roughness_ratio = np.array([0.05])
     convective_cooling.wind_speed = np.array([100.0])
     expected_result = np.array([115.5214])
 
-    result = convective_cooling._nu_forced(Tf, nu)
+    result = convective_cooling._nusselt_forced(Tf, nu)
 
     np.testing.assert_allclose(result, expected_result, rtol=1e-5)
 
@@ -107,12 +107,12 @@ def test_nu_forced_boundary_conditions(convective_cooling, expected_type):
     conv_cool_instances,
     ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"],
 )
-def test_nu_natural_single_value(convective_cooling, expected_type):
+def test_nusselt_natural_single_value(convective_cooling, expected_type):
     air_temperature = 50.0
     Td = 25.0
     nu = Air.kinematic_viscosity(air_temperature)
 
-    result = convective_cooling._nu_natural(air_temperature, Td, nu)
+    result = convective_cooling._nusselt_natural(air_temperature, Td, nu)
 
     assert isinstance(result, expected_type)
     np.testing.assert_allclose(result, [3.42404], rtol=1e-4)
@@ -123,12 +123,12 @@ def test_nu_natural_single_value(convective_cooling, expected_type):
     conv_cool_instances,
     ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"],
 )
-def test_nu_natural_array_values(convective_cooling, expected_type):
+def test_nusselt_natural_array_values(convective_cooling, expected_type):
     air_temperature = np.array([50.0, 60.0, 70.0])
     Td = np.array([25.0, 35.0, 45.0])
     nu = Air.kinematic_viscosity(air_temperature)
 
-    result = convective_cooling._nu_natural(air_temperature, Td, nu)
+    result = convective_cooling._nusselt_natural(air_temperature, Td, nu)
 
     assert isinstance(result, np.ndarray)
     assert result.shape == air_temperature.shape
@@ -140,12 +140,12 @@ def test_nu_natural_array_values(convective_cooling, expected_type):
     conv_cool_instances,
     ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"],
 )
-def test_nu_natural_edge_case(convective_cooling, expected_type):
+def test_nusselt_natural_edge_case(convective_cooling, expected_type):
     air_temperature = 0.0
     Td = 0.0
     nu = Air.kinematic_viscosity(air_temperature)
 
-    result = convective_cooling._nu_natural(air_temperature, Td, nu)
+    result = convective_cooling._nusselt_natural(air_temperature, Td, nu)
 
     assert isinstance(result, expected_type)
     assert result == 0.0
@@ -156,12 +156,12 @@ def test_nu_natural_edge_case(convective_cooling, expected_type):
     conv_cool_instances,
     ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"],
 )
-def test_nu_natural_high_values(convective_cooling, expected_type):
+def test_nusselt_natural_high_values(convective_cooling, expected_type):
     air_temperature = 1000.0
     Td = 500.0
     nu = Air.kinematic_viscosity(air_temperature)
 
-    result = convective_cooling._nu_natural(air_temperature, Td, nu)
+    result = convective_cooling._nusselt_natural(air_temperature, Td, nu)
 
     assert isinstance(result, expected_type)
     np.testing.assert_allclose(result, [2.18853], rtol=1e-4)
diff --git a/test/unit/solver/test_slv1t.py b/test/unit/solver/test_slv1t.py
index ccc3ce38..e4e5c03e 100644
--- a/test/unit/solver/test_slv1t.py
+++ b/test/unit/solver/test_slv1t.py
@@ -11,11 +11,12 @@
 
 from thermohl import power
 from thermohl.solver.entities import PowerType, VariableType
+from thermohl.power.rte.solar_heating import SolarHeating as RteSolarHeating
 from thermohl.solver.slv1t import Solver1T
 
 
 @pytest.fixture
-def solver():
+def solver_args():
     args = {
         "max_len": lambda: 1,
         "transit": np.array([0]),
@@ -25,8 +26,6 @@ def solver():
         "ambient_pressure": np.array([101325]),
         "relative_humidity": np.array([50]),
         "precipitation_rate": np.array([0]),
-        "linear_mass": 1.0,
-        "heat_capacity": 1.0,
         "datetime_utc": datetime(2025, 1, 1, 0, tzinfo=timezone.utc),
         "latitude": np.array([48.0]),
         "longitude": np.array([2.3]),
@@ -35,6 +34,8 @@ def solver():
         "solar_absorptivity": np.array([0.5]),
     }
     cable = {
+        "linear_mass": 1.0,
+        "heat_capacity": 1.0,
         "outer_diameter": np.array([3.186e-02]),
         "core_diameter": np.array([0.000]),
         "outer_area": np.array([6.004e-4]),
@@ -47,14 +48,26 @@ def solver():
         "emissivity": np.array([0.8]),
     }
     args.update(cable)
+    return args
+
+
+@pytest.fixture
+def solver(solver_args):
+    return create_solver(solver_args)
 
+
+def create_solver(solver_args):
     joule = power.rte.joule_heating.JouleHeating
     solar = power.rte.solar_heating.SolarHeating
     convective = power.rte.convective_cooling.ConvectiveCooling
     radiative = power.rte.radiative_cooling.RadiativeCooling
 
     solver = Solver1T(
-        dic=args, joule=joule, solar=solar, convective=convective, radiative=radiative
+        dic=solver_args,
+        joule=joule,
+        solar=solar,
+        convective=convective,
+        radiative=radiative,
     )
     return solver
 
@@ -218,3 +231,175 @@ def test_steady_intensity_custom_params(solver):
     for value in result.values():
         assert isinstance(value, np.ndarray)
     assert VariableType.TRANSIT.value in result
+
+
+def test_reduced_intensity_scalar_using_default_args() -> None:
+    args = {
+        "max_len": lambda: 1,
+        VariableType.TRANSIT.value: 20,
+        "ambient_temperature": 25,
+        "wind_speed": 0,
+        "wind_azimuth": 30.0,
+        "ambient_pressure": 101325,
+        "relative_humidity": 50,
+        "precipitation_rate": 0,
+        "month": 1,
+        "day": 1,
+        "hour": 0,
+        "latitude": 48.0,
+        "longitude": 2.3,
+        "altitude": 50.0,
+        "cable_azimuth": 0.0,
+        "solar_absorptivity": 0.5,
+    }
+    cable = {
+        "linear_mass": 1.0,
+        "heat_capacity": 1.0,
+        "outer_diameter": 3.186e-02,
+        "core_diameter": 0.000,
+        "outer_area": 6.004e-4,
+        "core_area": 0.000,
+        "magnetic_coeff": 1.000,
+        "magnetic_coeff_per_a": 0.000,
+        "temperature_coeff_linear": 3.600e-3,
+        "temperature_coeff_quadratic": 8.000e-7,
+        "linear_resistance_dc_20c": 5.540e-5,
+        "emissivity": 0.8,
+    }
+    args.update(cable)
+    solver = create_solver(args)
+
+    initial_wind_attack_angle = solver.convective_cooling.wind_attack_angle
+    initial_transit = args[VariableType.TRANSIT.value]
+
+    result = solver.reduced_intensity(
+        measured_temperature_difference=10.0,
+        measured_intensity=360.0,
+    )
+
+    assert isinstance(result, np.float64)
+
+    # Check that solver args and power term attributes have not been changed
+    assert np.isnan(solver.args.measured_global_radiation)
+    assert solver.args.ambient_temperature == args["ambient_temperature"]
+    assert solver.args.wind_speed == args["wind_speed"]
+    assert solver.args.transit == args["transit"]
+    assert np.isclose(solver.args.wind_attack_angle, np.deg2rad(30.0))
+
+    assert solver.convective_cooling.wind_attack_angle == initial_wind_attack_angle
+    assert solver.joule_heating.transit == initial_transit
+    assert isinstance(solver.solar_heating, RteSolarHeating)
+
+
+def test_reduced_intensity_scalar_providing_custom_args() -> None:
+    args = {
+        "max_len": lambda: 1,
+        VariableType.TRANSIT.value: 20,
+        "ambient_temperature": 25,
+        "wind_speed": 0,
+        "wind_azimuth": 30.0,
+        "ambient_pressure": 101325,
+        "relative_humidity": 50,
+        "precipitation_rate": 0,
+        "month": 1,
+        "day": 1,
+        "hour": 0,
+        "latitude": 48.0,
+        "longitude": 2.3,
+        "altitude": 50.0,
+        "cable_azimuth": 0.0,
+        "solar_absorptivity": 0.5,
+    }
+    cable = {
+        "linear_mass": 1.0,
+        "heat_capacity": 1.0,
+        "outer_diameter": 3.186e-02,
+        "core_diameter": 0.000,
+        "outer_area": 6.004e-4,
+        "core_area": 0.000,
+        "magnetic_coeff": 1.000,
+        "magnetic_coeff_per_a": 0.000,
+        "temperature_coeff_linear": 3.600e-3,
+        "temperature_coeff_quadratic": 8.000e-7,
+        "linear_resistance_dc_20c": 5.540e-5,
+        "emissivity": 0.8,
+    }
+    args.update(cable)
+    solver = create_solver(args)
+
+    initial_wind_attack_angle = solver.convective_cooling.wind_attack_angle
+    initial_transit = args[VariableType.TRANSIT.value]
+
+    result = solver.reduced_intensity(
+        measured_temperature_difference=10.0,
+        measured_intensity=360.0,
+        ambient_temperature=25.0,
+        wind_speed=4.0,
+        solar_irradiance=800.0,
+        max_conductor_temperature=120.0,
+    )
+
+    assert isinstance(result, np.float64)
+
+    # Check that solver args and power term attributes have not been changed
+    assert np.isnan(solver.args.measured_global_radiation)
+    assert solver.args.ambient_temperature == args["ambient_temperature"]
+    assert solver.args.wind_speed == args["wind_speed"]
+    assert solver.args.transit == args["transit"]
+    assert np.isclose(solver.args.wind_attack_angle, np.deg2rad(30.0))
+
+    assert solver.convective_cooling.wind_attack_angle == initial_wind_attack_angle
+    assert solver.joule_heating.transit == initial_transit
+    assert isinstance(solver.solar_heating, RteSolarHeating)
+
+
+def test_reduced_intensity_array() -> None:
+    args = {
+        "max_len": lambda: 2,
+        VariableType.TRANSIT.value: np.array([0, 0]),
+        "ambient_temperature": np.array([25, 25]),
+        "wind_speed": np.array([0, 0]),
+        "wind_azimuth": np.array([0, 0]),
+        "wind_attack_angle": np.array([45.0, 45.0]),
+        "ambient_pressure": np.array([101325, 101325]),
+        "relative_humidity": np.array([50, 50]),
+        "precipitation_rate": np.array([0, 0]),
+        "month": np.array([1, 1]),
+        "day": np.array([1, 1]),
+        "hour": np.array([0, 0]),
+        "latitude": np.array([48.0, 48.0]),
+        "longitude": np.array([2.3, 2.3]),
+        "altitude": np.array([50.0, 50.0]),
+        "cable_azimuth": np.array([0.0, 0.0]),
+        "solar_absorptivity": np.array([0.5, 0.5]),
+    }
+    cable = {
+        "linear_mass": np.array([1.0, 1.0]),
+        "heat_capacity": np.array([1.0, 1.0]),
+        "outer_diameter": np.array([3.186e-02, 3.186e-02]),
+        "core_diameter": np.array([0.000, 0.000]),
+        "outer_area": np.array([6.004e-4, 6.004e-4]),
+        "core_area": np.array([0.000, 0.000]),
+        "magnetic_coeff": np.array([1.000, 1.000]),
+        "magnetic_coeff_per_a": np.array([0.000, 0.000]),
+        "temperature_coeff_linear": np.array([3.600e-3, 3.600e-3]),
+        "temperature_coeff_quadratic": np.array([8.000e-7, 8.000e-7]),
+        "linear_resistance_dc_20c": np.array([5.540e-5, 5.540e-5]),
+        "emissivity": np.array([0.8, 0.8]),
+    }
+    args.update(cable)
+
+    solver = create_solver(args)
+
+    result = solver.reduced_intensity(
+        measured_temperature_difference=np.array([10.0, 10.0]),
+        measured_intensity=np.array([360.0, 360.0]),
+    )
+
+    assert isinstance(result, np.ndarray)
+    assert len(result) == 2
+    assert not np.isnan(result[0])
+    assert not np.isnan(result[1])
+
+    # Assert solver wind attack angle has not been changed
+    assert np.allclose(solver.args.wind_attack_angle, args["wind_attack_angle"])
diff --git a/test/unit/test_utils.py b/test/unit/test_utils.py
index 7b5db6fd..7d7ab06d 100644
--- a/test/unit/test_utils.py
+++ b/test/unit/test_utils.py
@@ -8,7 +8,7 @@
 import numpy as np
 import pytest
 
-from thermohl.utils import bisect_v, quasi_newton_2d
+from thermohl.utils import bisect_v, quasi_newton, quasi_newton_2d
 
 _nprs = 3141592654
 
@@ -102,7 +102,97 @@ def f(x):
         )
 
 
-#
+def test_quasi_newton_wrong_tol() -> None:
+    def f(x):
+        return x**2 - 2
+
+    with pytest.raises(ValueError):
+        quasi_newton(
+            f,
+            x0=1.0,
+            tol=-1.0e-6,  # Invalid tolerance
+        )
+
+
+def test_quasi_newton_wrong_maxiter() -> None:
+    def f(x):
+        return x**2 - 2
+
+    with pytest.raises(ValueError):
+        quasi_newton(
+            f,
+            x0=1.0,
+            maxiter=0,  # Invalid max iterations
+        )
+
+
+def test_quasi_newton_wrong_rtol() -> None:
+    def f(x):
+        return x**2 - 2
+
+    with pytest.raises(ValueError):
+        quasi_newton(
+            f,
+            x0=1.0,
+            rtol=-1.0e-6,  # Invalid relative tolerance
+        )
+
+
+def test_quasi_newton_scalar_increasing() -> None:
+    def f(x):
+        return np.log(x)
+
+    root = quasi_newton(
+        f,
+        x0=0.5,
+    )
+    assert np.isclose(root, 1.0)
+
+
+def test_quasi_newton_scalar_decreasing() -> None:
+    def f(x):
+        return -np.log(x)
+
+    root = quasi_newton(
+        f,
+        x0=0.5,
+    )
+    assert np.isclose(root, 1.0)
+
+
+def test_quasi_newton_scalar_no_convergence() -> None:
+    def f(x):
+        return x**2 + 2
+
+    with pytest.raises(RuntimeError):
+        quasi_newton(
+            f,
+            x0=1.0,
+        )
+
+
+def test_quasi_newton_array_no_convergence() -> None:
+    def f(x: np.ndarray) -> np.ndarray:
+        return x**2 + np.array([1, 2])
+
+    with pytest.raises(RuntimeError):
+        quasi_newton(
+            f,
+            x0=np.array([1.0, 1.0]),
+        )
+
+
+def test_quasi_newton_array_mixed() -> None:
+    def f(x: np.ndarray) -> np.ndarray:
+        return x**2 + np.array([-1, 2])
+
+    result = quasi_newton(
+        f,
+        x0=np.array([1.0, 1.0]),
+    )
+    assert np.isclose(result[0], 1.0)
+
+
 def test_quasi_newton_2d_convergence():
     np.random.seed(_nprs)
     size = 10
diff --git a/thermohl-docs/docs/api-reference/parameters.md b/thermohl-docs/docs/api-reference/parameters.md
index 3ccff0bc..311951e3 100644
--- a/thermohl-docs/docs/api-reference/parameters.md
+++ b/thermohl-docs/docs/api-reference/parameters.md
@@ -32,7 +32,8 @@ units, default values and in which set of power terms they are used.
 | hour      | 12            | N/A        | yes           | yes          | yes          | yes         | hour of the day (float in[0, 24[)                            |
 | ambient_temperature        | 15            | celsius    | yes           | yes          | yes          | yes         | ambient temperature                                          |
 | wind_speed        | 0             | linear_mass.s⁻¹      | yes           | yes          | yes          | yes         | wind speed                                                   |
-| wind_azimuth        | 90            | degree     | yes           | yes          | yes          | yes         | wind_azimuth (regarding north)                                 |
+| wind_azimuth        | 90            | degree     | yes           | yes          | yes          | yes         | wind azimuth (regarding north)                                 |
+| wind_attack_angle        | None            | radian     | yes           | yes          | yes          | yes         | (absolute) angle between the wind and the cable. Must be comprised between 0 and π / 2. If provided, this overrides the angle computed from wind_azimuth and cable_azimuth.                                  |
 | albedo        | 0.8           | N/A        | yes           | no           | no           | no          | albedo                                                       |
 | turbidity        | 0.1           | N/A        | no            | yes          | no           | no          | coefficient for air pollution from 0 (clean) to 1 (polluted) |
 | transit   | 100           | A          | yes           | yes          | yes          | yes         | transit intensity                                            |
diff --git a/thermohl-docs/docs/examples/ex_ist_reduction.ipynb b/thermohl-docs/docs/examples/ex_ist_reduction.ipynb
new file mode 100644
index 00000000..dc8ffb5c
--- /dev/null
+++ b/thermohl-docs/docs/examples/ex_ist_reduction.ipynb
@@ -0,0 +1,82 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8c937b14",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from datetime import datetime\n",
+    "\n",
+    "import numpy as np\n",
+    "\n",
+    "from thermohl import solver\n",
+    "from thermohl.solver import HeatEquationType\n",
+    "from thermohl.solver.entities import VariableType\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "81557e4b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "args = {\n",
+    "    \"max_len\": lambda: 1,\n",
+    "    VariableType.TRANSIT: np.array([0]),\n",
+    "    \"ambient_temperature\": np.array([25]),\n",
+    "    \"wind_speed\": np.array([0]),\n",
+    "    \"wind_azimuth\": np.array([0]),\n",
+    "    \"ambient_pressure\": np.array([101325]),\n",
+    "    \"relative_humidity\": np.array([50]),\n",
+    "    \"precipitation_rate\": np.array([0]),\n",
+    "    \"linear_mass\": 1.0,\n",
+    "    \"heat_capacity\": 1.0,\n",
+    "    \"datetime_utc\": datetime(2026, 1, 1, 0, 0, 0),\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0cf8da43",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "slvr = solver.ieee(\n",
+    "    args,\n",
+    "    heat_equation=HeatEquationType.ONE_TEMPERATURE,\n",
+    ")\n",
+    "\n",
+    "slvr.reduced_intensity(\n",
+    "    measured_temperature_difference=10.0,\n",
+    "    measured_intensity=360.0,\n",
+    "    max_conductor_temperature=100.0,\n",
+    ")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "thermohl (3.13.12)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.13.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/thermohl-docs/docs/user-guide.md b/thermohl-docs/docs/user-guide.md
index e3f02648..060eaee9 100644
--- a/thermohl-docs/docs/user-guide.md
+++ b/thermohl-docs/docs/user-guide.md
@@ -154,6 +154,30 @@ specific model with three temperatures for the conductor :
 * the average temperature;
 * the core temperature.
 
+## IST reduction (intensity limit reduction)
+
+For most computations, cable sleeves are ignored.
+However, if a sleeve is faulty, it can overheat. ThermoHL enables the user to
+compute a reduced intensity limit (called reduced IST) so that the interface
+between the sleeve and the cable doesn't exceed a given maximum temperature.
+Input data are
+- the measured temperature difference between the hotspot (at the interface between sleeve and cable)
+and the cable far away from the sleeve
+- the measured transit.
+
+This computation uses the single temperature model.
+
+Please note that following default values differ from the defaults for other
+computations:
+
+| Parameter | Default value |
+|----|----|
+| ambient_temperature | 30.0 °C |
+| wind_speed | 0.6 m/s |
+| measured_global_radiation | 600.0 W/m² |
+
+TODO: precision on wind_speed and wind_angle if needed
+
 ### Uncertainty in cable temperature computation
 
 For the conductor temperature computation with the three-temperatures "legacy" solver,
diff --git a/thermohl-docs/mkdocs.yml b/thermohl-docs/mkdocs.yml
index d82047de..0c677863 100644
--- a/thermohl-docs/mkdocs.yml
+++ b/thermohl-docs/mkdocs.yml
@@ -36,10 +36,11 @@ nav:
   - Home: index.md
   - Getting Started: getting-started.md
   - Examples:
-      - Example of steady temperature with 1t solver and IEEE model: examples/ex_steady_1t.ipynb
-      - Example of steady temperature with 3t solver and CIGRE model: examples/ex_steady_3t.ipynb
-      - Example of transient temperature with 3t solver and OLLA model: examples/ex_transient_3t.ipynb
-      - Example of steady temperature with 1t solver, all models and all parameters: examples/ex_all_params.ipynb
+      - Steady temperature with 1t solver and IEEE model: examples/ex_steady_1t.ipynb
+      - Steady temperature with 3t solver and CIGRE model: examples/ex_steady_3t.ipynb
+      - Transient temperature with 3t solver and OLLA model: examples/ex_transient_3t.ipynb
+      - Steady temperature with 1t solver, all models and all parameters: examples/ex_all_params.ipynb
+      - IST reduction: examples/ex_ist_reduction.ipynb
   - Parameters: api-reference/parameters.md
   - User Guide: user-guide.md
   - Developer Guide: