clear
clc
fileMeasurement1 = "E:\Studium_WS2023-2024\Fasttube\Tire Model\TireModel_self\Model_2\ft24_pacejka_tyre-main\ft24_pacejka_tyre\sorted_data\Hoosier 43164 20.5x7.0-13 R20, 8 inch rim\cornering.mat";
fileMeasurement2 = "E:\Studium_WS2023-2024\Fasttube\Tire Model\TireModel_self\Model_2\ft24_pacejka_tyre-main\ft24_pacejka_tyre\sorted_data\Hoosier 43164 20.5x7.0-13 R20, 8 inch rim\drivebrake.mat";

%
% cornering

data1 = load(fileMeasurement1);
n=1:length(data1.IA);
IA = round(data1.IA,0);
P = round(data1.P,0);%in kPa
FZ = data1.FZ; % in N
RL = data1.RL; % in cm
ia_0 = [];
ia_2 = [];
ia_4 = [];
for i=1:length(IA)
    if IA(i) == 0
        ia_0 = [ia_0, i];
    elseif IA(i) == 2
        ia_2 = [ia_2, i];
    elseif IA(i) == 4
        ia_4 = [ia_4, i];
    else
        continue
    end
end


p_50 = [];
p_70 = [];
p_80 = [];
p_90 = [];
for i=1:length(P)
    if P(i)< 60
        p_50 = [p_50, i];
    elseif P(i)< 75 && P(i) >= 60
        p_70 = [p_70, i];
    elseif P(i)< 90 && P(i) >= 75
        p_80 = [p_80, i];
    elseif P(i) < 100 && P(i) >= 90
        p_90 = [p_90, i];
    else
        continue
    end
end

ia = {ia_0, ia_2, ia_4};
p = {p_50, p_70, p_80, p_90};
for i=1:length(ia)
    for j=1:length(p)
        id{i,j} = intersect(ia{i}, p{j});
        [RL_temp, id_temp] = sort(RL(id{i,j}));% id_temp starts from 1

        FZ_temp = FZ(id{i,j});
        FZ_temp = FZ_temp(id_temp);

        [A, B] = createFit1(RL_temp, FZ_temp);
        slope{i,j} = A.p1;        
    end
end

%% drive brake
%

data2 = load(fileMeasurement2);
n=1:length(data2.IA);
IA = round(data2.IA,0);
P = round(data2.P,0); %in kPa
FZ = data2.FZ,2; % in N
RL = data2.RL,2; % in cm
ia_0 = [];
ia_2 = [];
ia_4 = [];
for i=1:length(IA)
    if IA(i) == 0
        ia_0 = [ia_0, i];
    elseif IA(i) == 2
        ia_2 = [ia_2, i];
    elseif IA(i) == 4
        ia_4 = [ia_4, i];
    else
        continue
    end
end

p_50 = [];
p_70 = [];
p_80 = [];
p_90 = [];
for i=1:length(P)
    if P(i)< 60
        p_50 = [p_50, i];
    elseif P(i)< 75 && P(i) >= 60
        p_70 = [p_70, i];
    elseif P(i)< 90 && P(i) >= 75
        p_80 = [p_80, i];
    elseif P(i) < 100 && P(i) >= 90
        p_90 = [p_90, i];
    else
       disp("P")
       disp(i)
    end
end
ia = {ia_0, ia_2, ia_4};
p = {p_50, p_70, p_80, p_90};
for i=1:length(ia)
    for j=1:length(p)
        id{i,j} = intersect(ia{i}, p{j});
        [RL_temp, id_temp] = sort(RL(id{i,j}));
        FZ_temp = FZ(id{i,j});
        FZ_temp = FZ_temp(id_temp);

        [A, B] = createFit1(RL_temp, FZ_temp);
        slope2{i,j} = A.p1;        
    end
end

%
%% Pairs

p_cell=[50,70,80,90];
ia_cell = [0,2,4];
disp("inflation pressure in Kpa")
disp(p_cell)
disp("inc Angle")
ia_cell

%
%% Results: Vertical stiffness Cornering

% row: incl angle
% column: inflation pressure
slope
figure()
for i=1:3
    stiffness = [slope{i,1}, slope{i,2}, slope{i,3}, slope{i,4}];
    plot(p_cell, stiffness)
    slope{i}
    hold on
end
legend("IA=0", "IA=2", "IA=4");
xlabel("inflation pressure in kPa")
ylabel("vertical stiffness in N/cm, cornering")

%% Results: Vertical stiffness DriveBrake

% row: incl angle
% column: inflation pressure
slope2
figure()
for i=1:3
    stiffness = [slope2{i,1}, slope2{i,2}, slope2{i,3}, slope2{i,4}];
    plot(p_cell, stiffness)
    hold on
end
legend("IA=0", "IA=2", "IA=4");
xlabel("inflation pressure in kPa")
ylabel("vertical stiffness in N/cm, drivebrake")

%% Conclusion:
% Ideally the relationship between vertical stiffness and inflation pressure should be linear and unidirectionally increasing.
%  reference:
% A Machine to Study Vertical Tire Stiffness and Damping Coefficient on JSTOR ;
% Measuring the tire vertical stiffness of a 255/40R17 ultra-high performance 200TW tire - Bridgestone Potenza RE-71R | Kelvin Tse (kktse.github.io) - Search (bing.com)
% Drive and Brake is more suitable to get the vertical stiffness, because it's more like the testing setup of tire vertical dynamics.