diff --git a/Rolling resistance/ComputeRollingResistance.m b/Rolling resistance/ComputeRollingResistance.m
new file mode 100644
index 0000000000000000000000000000000000000000..907225bfac4aa811e0f4842d11dfd766d5c81f40
--- /dev/null
+++ b/Rolling resistance/ComputeRollingResistance.m	
@@ -0,0 +1,59 @@
+%% ComputeRollingResistance
+% Use the original castor wheel on treadmill data to compute the rolling
+% resistance of the wheel when rolling on this treadmill. This parameter is
+% then used in the CastorModel as the rolling resistance coefficient (f_r).
+%   Date last edit: 2024-06-23
+%   Author: Matto Leeuwis
+
+% Get Data folder
+path = "../Data/Test8_31-5-2022/"; 
+
+% Get files for each trial where the maximal weight was used
+files = [
+    "0rad_7.4kg_8.txt";
+    "0.05rad_7.4kg_8.txt";
+    "0.1rad_7.4kg_8.txt";
+    "0.15rad_7.4kg_8.txt";
+    "0.2rad_7.4kg_8.txt"
+    ];
+
+% Initialize empty arrays
+f_r = [];
+f_r_std = [];
+phi_x = [0, 0.05, 0.1, 0.15, 0.2];
+
+for i = 1:length(files)
+    data = readtable([path + files(i)]);
+    
+    % Assume Fx is the rolling resistance force and Fz is the normal force
+    Fx = data.Fx_N_;  % Longitudinal force (direction of travel)
+    Fz = data.Fz_N_;  % Normal force
+    
+    % Calculate the Rolling Resistance Coefficient (f_r)
+    f_r(i) = -mean(Fx) / mean(Fz);
+
+    % Calculate standard deviation
+    f_r_i = -Fx./Fz;
+    f_r_std(i) = std(f_r_i);
+
+    % Display the rolling resistance coefficient
+    disp(['Calculated Rolling Resistance Coefficient (f_r): ', num2str(f_r(i))])
+end
+
+f_r;
+
+%% Perform linear regression
+mdl = fitlm(phi_x, abs(f_r))
+figure; plot(mdl)
+
+%% Print outcomes
+disp("Mean and std: ")
+for i = 1:length(f_r)
+    fprintf("%.2f rad:    %.3f, %.3f \n", phi_x(i), f_r(i), f_r_std(i))
+end
+
+fprintf("\nThe mean rolling resistance coefficient factor is %.4f\n", mean(f_r))
+
+
+
+