clear;clc;
delta = 2e-3;
F = 0:300;
b = 5e-3;
h = 20e-3;
L1 = 200e-3;
L2 = 250e-3;
kn = 300e3;
E = 210e9;
a = 2e-3;
ab = a/b;

syms F_p x1 x2 x
% kontrola
% double(solve(int(((F_p*x1)/(E*Jy(b,h))) * x1, x1, 0, L2) + int(( ((F_p*(x2+L2)) - (F*x2))/(E*Jy(b,h)) ) * (x2+L2), x2, 0, L1) == (-1/kn)*F_p - delta, F_p))

% vycisleni meho vypoctu
Mo = zeros([1 301]);
Fp = (((1/6)*((2*F*L1^3)+(3*F*L2*L1^2)))-delta*E*Jy(b,h)) / ((E*Jy(b,h))/(kn) + ((L1+L2)^3)/3)
for i = 1:length(F)
    if double(Fp(i)) >= 0
        Mo(i) = Fp(i)*(L1+L2) - F(i)*(L1);
    else
        Mo(i) = F(i) * L1;
    end 
end

Mo
sigma = (Mo * 0.5 * b)/(Jy(b,h))
y = 1.122 - 1.4*(ab) + 7.33*(ab^2) - 13.08*(ab^3) + 14*(ab^4)

plot(sigma * sqrt(pi * a) * y) 
xlim([0 300])
xlabel('F [N]')
ylabel('K_I [Pa?]')

xlim([0 300])




function J = Jy(b, h)
    J = ((b)*(h^3))/12;
end