diff --git a/src/fatpy/core/damage_cumulation/damage_cumulation_palmgren_meiner.py b/src/fatpy/core/damage_cumulation/damage_cumulation_palmgren_meiner.py
new file mode 100644
index 0000000..c2330d1
--- /dev/null
+++ b/src/fatpy/core/damage_cumulation/damage_cumulation_palmgren_meiner.py
@@ -0,0 +1,88 @@
+"""Damage Accumulation Rule according to Palmgren-Miner.
+
+Resources:
+    [1] Graphical interpretation of the change in the S-N curve based on the
+        chosen version can be found e.g. in this open-access paper:
+        http://dx.doi.org/10.5545/sv-jme.2013.1348
+    [2] Miner, M. A. (1945). Cumulative damage in fatigue. Journal of Applied
+        Mechanics, 12(3), 159-164.
+"""
+
+# import numpy as np
+# from numpy.typing import NDArray
+
+
+def damage_cumulation_elementary(
+    slope_k: float,
+    constant: float,
+    sig: float,
+    number_occurrences: int,
+) -> float:
+    r"""Elementary version of Palmgren-Miner linear damage accumulation.
+
+    The same slope k of the S-N curve below and above the fatigue limit.
+
+    ??? abstract "Math Equations"
+        $$
+        D = n/N = n\,\frac{\sigma^k}{C}
+        $$
+
+    """
+    total_occurrences: float = constant / sig**slope_k
+
+    damage: float = number_occurrences / total_occurrences
+
+    return damage
+
+
+def damage_cumulation_basic(
+    slope_k: float,
+    constant: float,
+    sig_fl: float,
+    sig: float,
+    number_occurrences: int,
+) -> float:
+    """Basic version of Palmgren-Miner linear damage accumulation.
+
+    The S-N curve gets horizontal at the fatigue limit, no damage for stresses beneath.
+    Otherwise elementary damage is calculated.
+    """
+    if sig < sig_fl:
+        damage = 0.0
+    else:
+        damage = damage_cumulation_elementary(
+            slope_k, constant, sig, number_occurrences
+        )
+
+    return damage
+
+
+def damage_cumulation_haibach(
+    slope_k: float,
+    constant: float,
+    sig_fl: float,
+    sig: float,
+    number_occurrences: int,
+) -> float:
+    r"""Haibach version of Palmgren-Miner linear damage accumulation.
+
+    the original slope_k is modified below fatigue limit to 2*slope_k-1.
+
+    ??? abstract "Math Equations"
+        $$
+        D = \frac{n}{C}\,\frac{\sigma^{2k-1}}{\sigma_\mathrm{FL}^{k-1}}
+        $$
+
+    """
+    if sig < sig_fl:
+        damage: float = (
+            number_occurrences
+            * sig ** (2 * slope_k - 1)
+            / (constant * sig_fl ** (slope_k - 1))
+        )
+    else:
+        damage = damage_cumulation_elementary(
+            slope_k, constant, sig, number_occurrences
+        )
+
+    return damage
diff --git a/src/fatpy/core/energy_life/damage_parameters.py b/src/fatpy/core/energy_life/damage_parameters.py
new file mode 100644
index 0000000..32ae385
--- /dev/null
+++ b/src/fatpy/core/energy_life/damage_parameters.py
@@ -0,0 +1,126 @@
+"""Damage parameters calculation methods for the energy-life."""
+
+import numpy as np
+from numpy.typing import NDArray
+from scipy.optimize import root_scalar
+from scipy.optimize import RootResults
+
+
+def calc_life_swt(
+    N: NDArray[np.float64],  # number of cycles
+    fat_strength_coef: NDArray[np.float64],
+    fat_strength_exp: NDArray[np.float64],
+    fat_ductility_coef: NDArray[np.float64],
+    fat_ductility_exp: NDArray[np.float64],
+    young_modulus: NDArray[np.float64],
+    p_swt: NDArray[np.float64],
+) -> NDArray[np.float64]:
+    """Function for root finding in SWT calculation."""
+    sol: NDArray[np.float64] = (
+        p_swt**2
+        - fat_strength_coef**2 * (2 * N) ** (2 * fat_strength_exp)
+        - young_modulus
+        * fat_ductility_coef
+        * fat_strength_coef
+        * (2 * N) ** (fat_strength_exp + fat_ductility_exp)
+    )
+    return sol
+
+
+def calc_dmg_param_swt(
+    elastic_modulus: NDArray[np.float64],
+    strain_amp: NDArray[np.float64],
+    mean_stress: NDArray[np.float64],
+    stress_amp: NDArray[np.float64],
+) -> NDArray[np.float64]:
+    """Calculate the SWT damage parameter."""
+    p_swt: NDArray[np.float64] = np.sqrt(
+        elastic_modulus * strain_amp * (mean_stress + stress_amp)
+    )
+    return p_swt
+
+
+# def swt(
+#     en_curve_parameters: dict[str, NDArray[np.float64]],
+#     stress_strain_values: dict[str, NDArray[np.float64]],
+#     N_0: float = 1.0,
+# ) -> int:
+#     """Calculate the number of cycles to failure according to SWT criterion."""
+#     n_values = swt_array(en_curve_parameters, stress_strain_values, N_0=N_0)
+#     return int(n_values.item())
+
+
+def swt(
+    en_curve_parameters: dict[str, NDArray[np.float64]],
+    stress_strain_values: dict[str, NDArray[np.float64]],
+    N_0: float = 1.0,
+) -> NDArray[np.int64]:
+    """Calculate SWT cycles to failure for array-like inputs.
+
+    All input values are broadcast to a common shape and solved elementwise.
+    """
+    fat_strength_coef = np.asarray(
+        en_curve_parameters["fat_strength_coef"], dtype=np.float64
+    )
+    fat_strength_exp = np.asarray(
+        en_curve_parameters["fat_strength_exp"], dtype=np.float64
+    )
+    fat_ductility_coef = np.asarray(
+        en_curve_parameters["fat_ductility_coef"], dtype=np.float64
+    )
+    fat_ductility_exp = np.asarray(
+        en_curve_parameters["fat_ductility_exp"], dtype=np.float64
+    )
+    elastic_modulus = np.asarray(
+        en_curve_parameters["elastic_modulus"], dtype=np.float64
+    )
+
+    strain_amp = np.asarray(stress_strain_values["strain_amp"], dtype=np.float64)
+    mean_stress = np.asarray(stress_strain_values["mean_stress"], dtype=np.float64)
+    stress_amp = np.asarray(stress_strain_values["stress_amp"], dtype=np.float64)
+
+    (
+        fat_strength_coef,
+        fat_strength_exp,
+        fat_ductility_coef,
+        fat_ductility_exp,
+        elastic_modulus,
+        strain_amp,
+        mean_stress,
+        stress_amp,
+    ) = np.broadcast_arrays(
+        fat_strength_coef,
+        fat_strength_exp,
+        fat_ductility_coef,
+        fat_ductility_exp,
+        elastic_modulus,
+        strain_amp,
+        mean_stress,
+        stress_amp,
+    )
+
+    if np.any(stress_amp <= np.abs(mean_stress)):
+        raise ValueError("SWT is only valid for stress_amp > |mean_stress|.")
+
+    p_swt = calc_dmg_param_swt(elastic_modulus, strain_amp, mean_stress, stress_amp)
+    roots = np.empty_like(p_swt, dtype=np.float64)
+
+    for idx in np.ndindex(p_swt.shape):
+        solution: RootResults = root_scalar(
+            calc_life_swt,
+            args=(
+                fat_strength_coef[idx],
+                fat_strength_exp[idx],
+                fat_ductility_coef[idx],
+                fat_ductility_exp[idx],
+                elastic_modulus[idx],
+                p_swt[idx],
+            ),
+            x0=N_0,
+            method="newton",
+        )
+        if not solution.converged:
+            raise ValueError(f"SWT calculation did not converge at index {idx}.")
+        roots[idx] = solution.root
+
+    return roots.astype(np.int64)
diff --git a/src/fatpy/core/energy_life/demage_parameters.py b/src/fatpy/core/energy_life/demage_parameters.py
deleted file mode 100644
index 03f895d..0000000
--- a/src/fatpy/core/energy_life/demage_parameters.py
+++ /dev/null
@@ -1 +0,0 @@
-"""Damage parameters calculation methods for the energy-life."""
diff --git a/tests/damage_cumulation/test_damage_cumulation.py b/tests/damage_cumulation/test_damage_cumulation.py
new file mode 100644
index 0000000..b40e421
--- /dev/null
+++ b/tests/damage_cumulation/test_damage_cumulation.py
@@ -0,0 +1,105 @@
+"""Test functions for damage accumulation rules."""
+
+import pytest
+import numpy as np
+# from numpy.typing import NDArray
+
+from fatpy.core.damage_cumulation import damage_cumulation_palmgren_meiner as dcpm
+
+
+@pytest.fixture
+def damage_cumulation_parameters() -> dict[str, float]:
+    """Fixture providing parameters for damage cumulation tests.
+
+    Returns:
+        dict[str, float]: Parameters including slope_k, constant, sig_fl.
+    """
+    params = {
+        "slope_k": 5.0,
+        "constant": 1e17,  # 1e15 is on the internet
+        "sig_fl": 137.97,
+    }
+    return params
+
+
+@pytest.fixture
+def fatigue_load_low() -> tuple[float, int]:
+    """Fixture providing a sample fatigue load.
+
+    Returns:
+        tuple[float, int]: Sample stress and number of occurrences.
+    """
+    return 150.0, 5000
+
+
+@pytest.fixture
+def fatigue_load_hi() -> tuple[float, int]:
+    """Fixture providing a sample fatigue load.
+
+    Returns:
+        tuple[float, int]: Sample stress and number of occurrences.
+    """
+    return 110.0, 100000
+
+
+def test_damage_cumulation_elementary(
+    damage_cumulation_parameters: dict[str, float],
+    fatigue_load_low: tuple[float, int],
+    fatigue_load_hi: tuple[float, int],
+) -> None:
+    """Elementary version
+    the same slope k of the S-N curve below and above the fatigue limit
+    """
+    slope_k = damage_cumulation_parameters["slope_k"]
+    constant = damage_cumulation_parameters["constant"]
+    sig_low, n_low = fatigue_load_low
+    sig_hi, n_hi = fatigue_load_hi
+
+    d_low = dcpm.damage_cumulation_elementary(slope_k, constant, sig_low, n_low)
+    d_hi = dcpm.damage_cumulation_elementary(slope_k, constant, sig_hi, n_hi)
+
+    assert np.around(d_low, decimals=4) == 0.0038
+    assert np.around(d_hi, decimals=4) == 0.0161
+
+
+def test_damage_cumulation_basic(
+    damage_cumulation_parameters: dict[str, float],
+    fatigue_load_low: tuple[float, int],
+    fatigue_load_hi: tuple[float, int],
+) -> None:
+    """Basic version
+    the S-N curve gets horizontal at the fatigue limit,
+    no damage for stresses beneath
+    """
+    slope_k = damage_cumulation_parameters["slope_k"]
+    constant = damage_cumulation_parameters["constant"]
+    sig_fl = damage_cumulation_parameters["sig_fl"]
+    sig_low, n_low = fatigue_load_low
+    sig_hi, n_hi = fatigue_load_hi
+
+    d_low = dcpm.damage_cumulation_basic(slope_k, constant, sig_fl, sig_low, n_low)
+    d_hi = dcpm.damage_cumulation_basic(slope_k, constant, sig_fl, sig_hi, n_hi)
+
+    assert np.around(d_low, decimals=4) == 0.0038
+    assert d_hi == 0.0
+
+
+def test_damage_cumulation_haibach(
+    damage_cumulation_parameters: dict[str, float],
+    fatigue_load_low: tuple[float, int],
+    fatigue_load_hi: tuple[float, int],
+) -> None:
+    """Haibach version
+    the original slope_k is modified below fatigue limit to 2*slope_k-1
+    """
+    slope_k = damage_cumulation_parameters["slope_k"]
+    constant = damage_cumulation_parameters["constant"]
+    sig_fl = damage_cumulation_parameters["sig_fl"]
+    sig_low, n_low = fatigue_load_low
+    sig_hi, n_hi = fatigue_load_hi
+
+    d_low = dcpm.damage_cumulation_haibach(slope_k, constant, sig_fl, sig_low, n_low)
+    d_hi = dcpm.damage_cumulation_haibach(slope_k, constant, sig_fl, sig_hi, n_hi)
+
+    assert np.around(d_low, decimals=4) == 0.0038
+    assert np.around(d_hi, decimals=5) == 0.00651
diff --git a/tests/energy_life/test_damage_parameters.py b/tests/energy_life/test_damage_parameters.py
new file mode 100644
index 0000000..2ef391c
--- /dev/null
+++ b/tests/energy_life/test_damage_parameters.py
@@ -0,0 +1,120 @@
+"""Test functions for damage parameters."""
+
+import pytest
+import numpy as np
+# from numpy.typing import NDArray
+
+from fatpy.core.energy_life import damage_parameters as dp
+
+
+@pytest.fixture
+def en_curve_parameters() -> dict[str, float]:
+    """Parameters of the e-N curve in the form of Manson-Coffin and Basquin equation.
+
+    Returns:
+        dict[str, float]: Parameters including:
+          fat_strength_coef: Manson-Coffin and Basquin equation fatigue
+                                strength coefficient
+          fat_ductility_coef: Manson-Coffin and Basquin equation fatigue
+                                ductility coefficient
+          fat_strength_exp: Manson-Coffin and Basquin equation fatigue
+                                strength exponent
+          fat_ductility_exp: Manson-Coffin and Basquin equation fatigue
+                                ductility exponent
+          elastic_modulus: Young's / Elastic modulus
+    """
+    params = {
+        "fat_strength_coef": 475.4,
+        "fat_ductility_coef": 0.612,
+        "fat_strength_exp": -0.078,
+        "fat_ductility_exp": -0.62,
+        "elastic_modulus": 162000.0,
+    }
+    return params
+
+
+@pytest.fixture
+def stress_strain_values() -> dict[str, float]:
+    """Stress / Strain values.
+
+    Returns:
+        dict[str, float]: Parameters including:
+          strain_amp: Strain amplitude
+          stress_amp: Stress amplitude
+          mean_stress: Mean stress
+    """
+    params = {"strain_amp": 0.0135, "stress_amp": 290.0, "mean_stress": 10.0}
+    return params
+
+
+@pytest.fixture
+def stress_strain_values_array() -> dict[str, np.ndarray]:
+    """Stress / Strain values as arrays for testing swt with array inputs.
+
+    Returns:
+        dict[str, np.ndarray]: Parameters including:
+          strain_amp: Strain amplitude
+          stress_amp: Stress amplitude
+          mean_stress: Mean stress
+    """
+    params = {
+        "strain_amp": np.array([0.0135, 0.0135]),
+        "stress_amp": np.array([290.0, 290.0]),
+        "mean_stress": np.array([10.0, 10.0]),
+    }
+    return params
+
+
+@pytest.fixture
+def invalid_stress_condition_array() -> dict[str, np.ndarray]:
+    """Stress / Strain values that violate SWT validity condition
+    (stress_amp <= abs(mean_stress)).
+    """
+
+    params = {
+        "strain_amp": np.array([0.0135, 0.0135]),
+        "stress_amp": np.array([10.0, 20.0]),
+        "mean_stress": np.array([10.0, 10.0]),
+    }
+    return params
+
+
+def test_swt(
+    en_curve_parameters: dict[str, float], stress_strain_values: dict[str, float]
+) -> None:
+    """Tests Smith-Watson-Topper (SWT) damage parameter for mean stress correction
+    in strain-life.
+
+    """
+
+    n = dp.swt(en_curve_parameters, stress_strain_values)
+
+    elastic_modulus: float = en_curve_parameters["elastic_modulus"]
+    strain_amp: float = stress_strain_values["strain_amp"]
+    mean_stress: float = stress_strain_values["mean_stress"]
+    stress_amp: float = stress_strain_values["stress_amp"]
+    p_swt = dp.calc_dmg_param_swt(elastic_modulus, strain_amp, mean_stress, stress_amp)
+
+    assert n == 278
+    assert p_swt == 810.0
+
+
+def test_swt_array_returns_elementwise_cycles(
+    en_curve_parameters: dict[str, float],
+    stress_strain_values_array: dict[str, np.ndarray],
+) -> None:
+    """swt_array should solve SWT for each broadcasted input entry."""
+
+    n_values = dp.swt(en_curve_parameters, stress_strain_values_array)
+
+    assert np.array_equal(n_values, np.array([278, 278]))
+
+
+def test_swt_array_raises_for_invalid_stress_condition(
+    en_curve_parameters: dict[str, float],
+    invalid_stress_condition_array: dict[str, np.ndarray],
+) -> None:
+    """SWT is not valid when stress_amp <= abs(mean_stress)."""
+
+    with pytest.raises(ValueError, match="SWT is only valid"):
+        dp.swt(en_curve_parameters, invalid_stress_condition_array)