L=27.9;
B=7.2;
T=2.65;
H=3.49;
delta=0.55;
alpha=0.80;
beta=0.84;
kz=0.6;
h0=2*B*sqrt(((1.017+0.023)*alpha/(alpha+delta))*((1.06+0.05)*alpha^2/12/delta))-kz*H;
te=[0 10 20 30 40 50 60 70 80 90];
Hte1=[0 0.050 0.387 0.840 1.279 1.365 1.056 0.583 0.210 0];
Hte2=[0 -0.036 -0.241 -0.556 -0.722 -0.513 0.026 0.603 0.935 1.000];
Hte3=[0 0.151 0.184 0.081 -0.069 -0.155 -0.135 -0.062 -0.010 0];
Hte4=[0 0.010 0.062 0.135 0.155 0.069 -0.081 -0.184 -0.151 0];
for i=1:90;
fte1(i)=interp1(te,Hte1,i);
fte2(i)=interp1(te,Hte2,i);
fte3(i)=interp1(te,Hte3,i);
fte4(i)=interp1(te,Hte4,i);
lst(i)=0.5*B*(1-0.972*T/H)*fte1(i)+0.64*(1-1.032*T/H)*H*fte2(i)+1/11.4*...
(alpha*B)^2/delta/T*fte3(i)+1/11.4*(alpha*B)^2/delta/T*((0.64*(1-1.032*T/H)*H)/...
(0.5*B*(1-0.972*T/H)))^3*fte4(i)-(kz*H-alpha/(alpha+delta)*T)*sin(i*pi/180);
end
ls=lst(max(lst)==lst);
tkr=find((lst)==ls);
for i=1:120;
fte1(i)=interp1(te,Hte1,i);
fte2(i)=interp1(te,Hte2,i);
fte3(i)=interp1(te,Hte3,i);
fte4(i)=interp1(te,Hte4,i);
lst1(i)=0.5*B*(1-0.972*T/H)*fte1(i)+0.64*(1-1.032*T/H)*H*fte2(i)+1/11.4*...
(alpha*B)^2/delta/T*fte3(i)+1/11.4*(alpha*B)^2/delta/T*((0.64*(1-1.032*T/H)*H)/...
(0.5*B*(1-0.972*T/H)))^3*fte4(i)-(kz*H-alpha/(alpha+delta)*T)*sin(i*pi/180);
if lst1(i)<0
disp(['u1=', num2str(lst1(i))]);
u2=find((lst1)==(lst1(i)));
break
end
end
for i=10:90;
fte1(i)=interp1(te,Hte1,i);
fte2(i)=interp1(te,Hte2,i);
fte3(i)=interp1(te,Hte3,i);
fte4(i)=interp1(te,Hte4,i);
lstd(i)=0.5*B*(1-0.972*T/H)*fte1(i)+0.64*(1-1.032*T/H)*H*fte2(i)+1/11.4*...
(alpha*B)^2/delta/T*fte3(i)+1/11.4*(alpha*B)^2/delta/T*((0.64*(1-1.032*T/H)*H)/...
(0.5*B*(1-0.972*T/H)))^3*fte4(i)-(kz*H-alpha/(alpha+delta)*T)*sin(i*pi/180);
end
ls0=0;
ls10=lstd(10);
ls20=lstd(20);
ls30=lstd(30);
ls40=lstd(40);
ls50=lstd(50);
ls60=lstd(60);
ls70=lstd(70);
ls80=lstd(80);
ls90=lstd(90);
X=0.0873;
l0=ls0*X;
l10=ls10*X;
l20=(2*ls10+ls20)*X;
l30=(2*ls10+2*ls20+ls30)*X;
l40=(2*ls10+2*ls20+2*ls30+ls40)*X;
l50=(2*ls10+2*ls20+2*ls30+2*ls40+ls50)*X;
l60=(2*ls10+2*ls20+2*ls30+2*ls40+2*ls50+ls60)*X;
l70=(2*ls10+2*ls20+2*ls30+2*ls40+2*ls50+2*ls60+ls70)*X;
l80=(2*ls10+2*ls20+2*ls30+2*ls40+2*ls50+2*ls60+2*ls70+ls80)*X;
l90=(2*ls10+2*ls20+2*ls30+2*ls40+2*ls50+2*ls60+2*ls70+2*ls80+ls90)*X;
kr=[0 10 20 30 40 50 60 70 80 90];
M1=[ls0 ls10 ls20 ls30 ls40 ls50 ls60 ls70 ls80 ls90];
M2=[l0 l10 l20 l30 l40 l50 l60 l70 l80 l90];
Fp0 = griddedInterpolant(kr,M1);
funp0 = @(t) Fp0(t);
Fp3=griddedInterpolant(kr,M2);
funp3 = @(t) Fp3(t);
plot(kr, funp0(kr),kr,funp3(kr),'linewidth',1.5);
grid on
xlabel('heel angel')
ylabel('arm')
legend('ls','ld','Location','northwest')
title('Diagram stability')
hold on;
