clc
clear all
close all

% This code is used to plot the ISO 8608 input road profile Classes A to E
% for both non-suppressed dataset and harmonic-suppressed dataset as well
% as the associated body and tyre displacements.

load suspension.mat
roadx = {hxA, hxB, hxC, hxD, hxE};
roads = {hsA, hsB, hsC, hsD, hsE};

bodyx = {bodyxsuppressedA, bodyxsuppressedB, bodyxsuppressedC, bodyxsuppressedD, bodyxsuppressedE};
bodys = {bodysuppressedA, bodysuppressedB, bodysuppressedC, bodysuppressedD, bodysuppressedE};

tyrex = {tirexsuppressedA, tirexsuppressedB, tirexsuppressedC, tirexsuppressedD, tirexsuppressedE};
tyres = {tiresuppressedA, tiresuppressedB, tiresuppressedC, tiresuppressedD, tiresuppressedE};

roadclass = {'A','B','C','D','E'};
linecolour = {'#0066cc','#cc6600','#ff9933','#663399','#669933'};


L = 2000;
B = 0.02;
x = 0:B:L-B;
N = 100000;
dn = 1/L;


for i = 1:5
    % No harmonic suppressed
    figure(i)
    plot(x,roadx{i},'Color',linecolour{i},'LineWidth',1)
    title(['Road profile elevation Class ' roadclass{i}])
    xlabel('Horizontal Distance [m]')
    ylabel('Vertical Displacement [m]')

    figure(222)    
    loglog([dn:dn:dn*N],(abs(fft(roadx{i})).^2)/N,'.','LineWidth',2)
    hold on
    grid
    title('Displacement PSD')
    xlabel('Spatial Frequency [cycles/m]')
    ylabel('Road profile elevation PSD [m^3]')
    xlim([1e-2 10])
    ylim([1e-10 1])
    legend('A','B','C','D','E','Location','Best')
    
    % With harmonic suppressed

    figure(i+10)
    plot(x,roads{i},'Color',linecolour{i},'LineWidth',1)
    title(['Road profile elevation Class ' roadclass{i}])
    xlabel('Horizontal Distance [m]')
    ylabel('Vertical Displacement [m]')
 
    figure(444)    
    loglog([dn:dn:dn*N],(abs(fft(roads{i})).^2)/N,'.','LineWidth',2)
    hold on
    grid
    title('Displacement PSD')
    xlabel('Spatial Frequency [cycles/m]')
    ylabel('Road profile elevation PSD [m^3]')
    xlim([1e-2 10])
    ylim([1e-10 1])
    legend('A','B','C','D','E','Location','Best')
end


%% Plotting output body and tyre displacements
% close all
Ts = 0.001;
N = 100000;
Nt = 200000;
t = 0:Ts:(100000-1)*Ts;
bodyxcut = {bodyxsuppressedA(N+1:Nt), bodyxsuppressedB(N+1:Nt), bodyxsuppressedC(N+1:Nt),...
    bodyxsuppressedD(N+1:Nt), bodyxsuppressedE(N+1:Nt)};
tyrexcut = {tirexsuppressedA(N+1:Nt), tirexsuppressedB(N+1:Nt), tirexsuppressedC(N+1:Nt),...
    tirexsuppressedD(N+1:Nt), tirexsuppressedE(N+1:Nt)};

bodyscut = {bodysuppressedA(N+1:Nt), bodysuppressedB(N+1:Nt), bodysuppressedC(N+1:Nt),...
    bodysuppressedD(N+1:Nt), bodysuppressedE(N+1:Nt)};
tyrescut = {tiresuppressedA(N+1:Nt), tiresuppressedB(N+1:Nt), tiresuppressedC(N+1:Nt),...
    tiresuppressedD(N+1:Nt), tiresuppressedE(N+1:Nt)};


for j = 1:5
    % No Harmonic Suppressed
    figure(j)
    plot(t,bodyxcut{j},'Color',linecolour{j},'LineWidth',1)
    xlabel('Time [s]')
    title(['Class ' roadclass{j}])
    ylabel('Body Displacement [m]')
    
    figure(j+10)
    plot(t,tyrexcut{j},'Color',linecolour{j},'LineWidth',1)
    xlabel('Time [s]')
    title(['Class ' roadclass{j}])
    ylabel('Tyre Displacement [m]')
    
    figure(666)
    loglog([dn:dn:dn*N],(abs(fft(bodyx{j}(1:N))).^2)/N,'LineWidth',2)
    hold on
    grid    
    xlabel('Spatial Frequency, cycles/m')
    ylabel('Body Displacement PSD [m^3]')
    xlim([1e-2 10])

    figure(777)
    loglog([dn:dn:dn*N],(abs(fft(tyrex{j}(1:N))).^2)/N,'LineWidth',2)
    hold on
    grid
    xlabel('Spatial Frequency [cycles/m]')
    ylabel('Tyre Displacement PSD [m^3]')
    xlim([1e-2 10])
 
%     Suppressed Harmonics
    figure(j+20)
    plot(t,bodyscut{j},'Color',linecolour{j},'LineWidth',1)    
    title(['Class ' roadclass{j}])
    ylabel('Displacement [m]')
    
    figure(j+30)
    plot(t,tyrescut{j},'Color',linecolour{j},'LineWidth',1)
    title(['Class ' roadclass{j}])
    ylabel('Displacement [m]')
 
    figure(888)
    loglog([dn:dn:dn*N],(abs(fft(bodys{j}(1:N))).^2)/N,'.','LineWidth',2)
    hold on
    grid    
    xlabel('Spatial Frequency [cycles/m]')
    ylabel('Road profile elevation PSD [m^3]')
    xlim([1e-2 10])
    legend('Class A','Class B','Class C','Class D','Class E','Location','Best')

    figure(999)
    loglog([dn:dn:dn*N],(abs(fft(tyres{j}(1:N))).^2)/N,'.','LineWidth',2)
    hold on
    grid
    xlabel('Spatial Frequency [cycles/m]')
    ylabel('Tyre Displacement PSD [m^3]')
    xlim([1e-2 10])
    legend('Class A','Class B','Class C','Class D','Class E','Location','Best')

end

