clear all; clc;
%Parametry
r=25;       %poloměr zvedacích válců v mm
h=40;       %výška uchycení válců v mm
m=2325;     %pohotovostní hmotnost vozu v kg
S_np=265;   %nominální šířka předních pneumatik v mm
S_nz=305;   %nominální šířka zadních pneumatik v mm
AR_p=35;    %poměrná výška profilu řpedních pneumatik
AR_z=30;    %poměrná výška profilu řpedních pneumatik
D_rp=21;    %průměr ráfků předních kol v palcích
D_rz=21;    %průměr ráfků předních kol v palcích
g=9.81;     %gravitační zrychlení
p_p=285;    %tlak předních pneumatik v kPa
p_z=305;    %tlak zadních pneumatik v kPa

G=m*g/4;

D_p=S_np*AR_p/50+D_rp*25.4;
D_z=S_nz*AR_z/50+D_rz*25.4;
R_p=D_p/2;
R_z=D_z/2;

k_p=(0.00028*p_p*sqrt((-0.004*AR_p+1.03)*S_np*D_p)+3.45)*g; 
k_z=(0.00028*p_z*sqrt((-0.004*AR_z+1.03)*S_nz*D_z)+3.45)*g;   

%Pneumatika se dotýká podložky
beta0_p=asin((R_p-h-G/k_p)/(R_p+r));
beta0_z=asin((R_z-h-G/k_z)/(R_z+r));
beta_p=linspace(beta0_p, pi/2, 5000);
beta_z=linspace(beta0_z, pi/2, 5000);

a_p=2/3*R_p-1./(3.*sin(beta_p)).*(G/k_p-R_p+h+r.*sin(beta_p));
a_z=2/3*R_z-1./(3.*sin(beta_z)).*(G/k_z-R_z+h+r.*sin(beta_z));

ksi1_p=R_p-h-(a_p+r).*sin(beta_p);
ksi1_z=R_z-h-(a_z+r).*sin(beta_z);
y_p=-ksi1_p;
y_z=-ksi1_z;

%Pneumatika se nedotýká podložky
for i=1:length(beta_p);
    if ksi1_p(i)<0
        a_p(i)=R_p-G./(2*k_p.*sin(beta_p(i)));
        y_p(i)=h+r.*sin(beta_p(i))+a_p(i)*sin(beta_p(i))-R_p;
        ksi1_p(i)=0;
    end
end

for i=1:length(beta_z);
    if ksi1_z(i)<0
        a_z(i)=R_z-G./(2*k_z.*sin(beta_z(i)));
        y_z(i)=h+r.*sin(beta_z(i))+a_z(i)*sin(beta_z(i))-R_z;
        ksi1_z(i)=0;
    end
end

F_p=k_p.*(R_p-a_p).*cos(beta_p);
F_z=k_z.*(R_z-a_z).*cos(beta_z);
L_p=2.*(a_p+r).*cos(beta_p);
L_z=2.*(a_z+r).*cos(beta_z);

figure('Name','Průběhy stahovací síly a zdvihu kola při uvažování poddajné pneumatiky','NumberTitle','off');
yyaxis left
plot(L_p,F_p, 'Linewidth', 1.5);
title({'Průběhy stahovací síly a zdvihu kola' 'při uvažování poddajné pneumatiky'});
xlabel('Rozteč zvedacích válců [mm]'); ylabel('Stahovací síla [N]');
hold on;
plot(L_z,F_z, 'Linewidth', 1.5);

yyaxis right
plot(L_p,y_p, 'Linewidth', 1.5);
ylabel ('Zdvih pneumatiky [mm]');
hold on;
plot(L_z,y_z, 'Linewidth', 1.5);

set(gca,'FontName','Technika Book','FontSize', 12,'XGrid','on','YGrid','on');
print -djpeg PrubehPoddajna

F_pmax=max(F_p)
F_zmax=max(F_z)
L_maxp=max(L_p)
L_maxz=max(L_z)