% BEAM THEORY
clear all; close all; clc;

% General properties

loadcase = 9;

L = 237e-3; % [m]
a = 158e-3; % [m]
b = L - a; % [m]
E = 69e9; % [Pa]
F = 49; % [N]
M = 10e-3; % [Nm]
q = F/L; % [N/m]

dx = 1e-3; % [m]
x = 0:dx:L; % [m]

% cross section area / area moment of inertia [m^4] / extreme fiber
% distance [m]
t = 15e-3; % thickness in [m]
B = 60e-3; % width in [m]


I = 1/12*B*t^3; %area moment of inertia [m^4]
u = t/2; %extreme fiber distance [m]
A=B*t;


%% Load case 1
if loadcase == 1;
C = 1;

dz = -(F.*x.^2).*(3.*L-x)./(6.*E.*I);
D  = ones(1,length(x)).*F;
MR = F.*(x-L);
sigma  = abs(F.*u.*(x-L)./I);

dz_max = -F*L^3/3/E/I;
xdz = L;
MR_max = -F*L;
xMR = 0;
sigma_max = abs(F*L*u/I);
xsigma = 0;
Cz = abs(-3*E*I/L^3);
xC = L;
end

%% Load case 2
if loadcase == 2;
C = 3;

dz = (-M.*x.^2)./(2.*E.*I);
D  = zeros(1,length(x));
MR = ones(1,length(x)).*M;
sigma  = abs(M.*u./I.*ones(1,length(x)));

dz_max = -M*L^2/2/E/I;
xdz = L;
MR_max = M;
xMR = -1;
sigma_max = abs(M*u/I);
xsigma = -1;
Cz = abs(-E*I/L);
xC = L;
end

%% Load case 3
if loadcase == 3;
C = 2;

dz = -(q.*x.^2./(24.*E.*I)).*(6.*L.^2-4.*L.*x+x.^2);
D  = q.*(L-x);
MR = -q.*(L-x).^2./2;
sigma  = abs(-q.*u.*(L-x).^2./2./I);

dz_max = -q*L^4/8/E/I;
xdz = L;
MR_max = -q*L^2/2;
xMR = 0;
sigma_max = abs(-q*L^2*u/2*I);
xsigma = 0;
Cz = abs(-8*E*I/L^4);
xC = L;
end

%% Load case 4
if loadcase == 4;
x1 = 0:dx:a;
x2 = a+dx:dx:L;
C = 1;
    
dz1 = -F.*a.*b.^2./(6.*E.*I).*((1+L./b).*(x1./L)-x1.^3./(a.*b.*L));
dz2 = -F.*a.^2.*b./(6.*E.*I).*((1+L./a).*(L-x2)./L-(L-x2).^3./(a.*b.*L));
D1 = ones(1,length(x1)).*F.*b./L;
D2  = ones(1,length(x2)).*(F.*b./L)-F;
MR1 = F.*b.*x1./(L);
MR2 = (F.*b.*x2./(L))-(F.*(x2-a));
sigma1 = F.*b.*x1.*u./(I.*L);
sigma2 = ((F.*b.*x2./(L))-(F.*(x2-a)))./I.*u;

dz = [dz1 dz2];
D = [D1 D2];
MR = [MR1 MR2];
sigma = abs([sigma1 sigma2]);

if a>=b
    dz_max = -F*b*(L^2-b^2)^1.5/(9*3^0.5*E*I*L);
    xdz = ((L^2-b^2)/3)^0.5;
elseif a<b
        dz_max = -F*a*(L^2-a^2)^1.5/(9*3^0.5*E*I*L);
        xdz = L-(((L^2-a^2)/3)^0.5);
end
MR_max = F*b*a/(a+b);
xMR = a;
sigma_max = abs(F*b*a*u/(I*(a+b)));
xsigma = a;
Cz = abs(-3*E*I*L/(a^2*b^2));
xC = a;
end

%% Load case 5
if loadcase == 5;
C = 3;

dz = -(M.*L.^2./(6.*E.*I)).*(x./L-(x./L).^3);
D  = ones(1,length(x)).*M./L;
MR = M./L.*x;
sigma  = abs(M.*x.*u./L/I);

dz_max = -M*L^2/(9*3^0.5*E*I);
xdz = L/3^0.5;
MR_max = M;
xMR = L;
sigma_max = abs(M*u/I);
xsigma = L;
Cz = abs(-3*E*I/L);
xC = L;
end

%% Load case 6
if loadcase == 6;
C = 2;

dz = -q.*L.^4./(24.*E.*I).*(x./L-2.*(x./L).^3+(x./L).^4);
D  = q.*L./2-q.*x;
MR = q.*x.*(L-x)/2;
sigma  = abs(q.*x.*u.*(L-x)./(2.*I));

dz_max = -5*q*L^4/(384*E*I);
xdz = L/2;
MR_max = q*L^2/8;
xMR = L/2;
sigma_max = abs(q*L^2*u/(8*I));
xsigma = L/2;
Cz = abs(-5*q*L^4/(384*E*I));
xC = L/2;
end

%% Load case 7
if loadcase == 7;
x1 = 0:dx:a;
x2 = a+dx:dx:L;
C = 1;
    
dz1 = F.*b.*a.^2/(6.*E.*I).*(x1./a-(x1./a).^3);
dz2 = -F.*((-x2)+a)./(6.*E.*I).*(a.*b-3.*b.*x2+x2.^2-2.*a.*x2+a.^2);
D1 = ones(1,length(x1)).*-(F.*b./a);
D2  = ones(1,length(x2)).*F;
MR1 = -(F.*b.*x1./a);
MR2 = -F.*((-x2)+a+b);
sigma1 = -F.*b.*x1.*u/I/a;
sigma2 = -(F.*((-x2)+a+b))./I.*u;

dz = [dz1 dz2];
D = [D1 D2];
MR = [MR1 MR2];
sigma = abs([sigma1 sigma2]);

dz_max1 = F*b*a^2/(9*3^0.5*E*I);
dz_max2 = -F*b^2*L/(3*E*I);
if abs(dz_max1)<abs(dz_max2)
    dz_max = dz_max2;
    xdz = L;
elseif abs(dz_max1)>=abs(dz_max2)
    dz_max = dz_max1;
    xdz = a/(3^0.5);
end
MR_max = -F*a;
xMR = a;
sigma_max = abs(-F*a*u/I);
xsigma = a;
Cz = abs(-3*E*I/b^2/L);
xC = L;
end

%% Load case 8
if loadcase == 8;
x1 = 0:dx:a;
x2 = a+dx:dx:L;
C = 4;
    
dz1 = q.*b.^2.*a.^2./(12.*E.*I).*(x1./a-(x1./a).^3);
dz2 = -q.*b.^4./(24.*E.*I).*(4.*(a./b).*((x2-a)./b)+6.*((x2-a)./b).^2-4.*((x2-a)./b).^3+((x2-a)./b).^4);
D1 = ones(1,length(x1)).*-q.*b.^2./(2.*a);
D2  = q.*b-q.*(x2-a);
MR1 = -q.*x1.*b.^2./(2.*a);
MR2 = -q.*(b.^2-2.*b.*x2+2.*b.*a+x2.^2-2.*x2.*a+a.^2)./2;
sigma1 = q.*x1.*u.*b.^2./(2.*a.*I);
sigma2 = q.*u.*(b.^2-2.*b.*x2+2.*b.*a+x2.^2-2.*x2.*a+a.^2)./(2.*I);

dz = [dz1 dz2];
D = [D1 D2];
MR = [MR1 MR2];
sigma = abs([sigma1 sigma2]);

dz_max1 = q*b^2*a^2/(18*3^0.5*E*I);
dz_max2 = -q*b^3*(4*a+3*b)/(24*E*I);
if abs(dz_max1)<abs(dz_max2)
    dz_max = dz_max2;
    xdz = L;
elseif abs(dz_max1)>=abs(dz_max2)
    dz_max = dz_max1;
    xdz = a/(3^0.5);
end
MR_max = -q*b^2/2;
xMR = a;
sigma_max = abs(-q*b^2*u/2/I);
xsigma = a;
Cz1 = abs(18*3^0.5*E*I/(b^2*a^2));
xC1 = a/(3^0.5);
Cz2 = abs(-24*E*I/(b^3*(4*a+3*b)));
xC2 = L;
end

%% Load case 9
if loadcase == 9;
x1 = 0:dx:a;
x2 = a+dx:dx:L;
C = 1;
    
dz1 = -F.*L.*b.^2./(4.*E.*I).*(a.*x1./L.^2-2./3.*(1+a./(2.*L)).*(x1./L).^3);
dz2 = F.*a.*(L-x2).^2./(12.*L.^3).*((2.*L.*a.^2-3.*L.^2.*x2+x2.*a.^2)./E./I);
D1 = ones(1,length(x1)).*F.*(b./L).^2.*(1+a./2./L);
D2  = ones(1,length(x2)).*-F.*(a./L).^2.*(1+b./2./L+3.*b./2./a);
MR1 = F.*x1.*(b./L).^2.*(1+a./2./L);
MR2 = F.*x2.*(b./L).^2.*(1+a./2./L)-F.*(x2-a);
sigma1 = F.*x1.*b.^2.*u.*(2.*L+a)./(2.*L.^3.*I);
sigma2 = F.*u.*(2.*x2.*L.*b.^2+b.^2.*x2.*a-2.*L.^3.*x2+2.*L.^3.*a)./(2.*L.^3.*I);

dz = [dz1 dz2];
D = [D1 D2];
MR = [MR1 MR2];
sigma = abs([sigma1 sigma2]);

x_max1 = L*((2*a*L+b*L-b*a)/(2*a*L+3*b*L+b*a));
x_max2 = L*((a/(2*L))/(1+a/(2*L)))^0.5;

if a<0.414*L
    dz_max = F.*a.*(L-x_max1).^2./(12.*L.^3).*((2.*L.*a.^2-3.*L.^2.*x_max1+x_max1.*a.^2)./E./I);
    xdz = L*((2*a*L+b*L-b*a)/(2*a*L+3*b*L+b*a));
    MR_max = F*a*b^2/L^2*(1+a/2/L);
    xMR = a;
    sigma_max = F*a*b^2*u*(1+a/2/L)/L^2/I;
    xsigma = a;
elseif a>=0.414*L
    dz_max = -F.*L.*b.^2./(4.*E.*I).*(a.*x_max2./L.^2-2./3.*(1+a./(2.*L)).*(x_max2./L).^3);
    xdz = L*((a/(2*L))/(1+a/(2*L)))^0.5;
    MR_max = -F*a*b/L*(1-b/2/L);
    xMR = L;
    sigma_max = abs(F*a*b*u*(1-b/2/L)/L/I);
    xsigma = L;
end
Cz = abs(12*E*I*L^3/(b^2*a^2*(3*L^2-2*a*L-a^2)));
xC = a;

end

%% Load case 10
if loadcase == 10;
x1 = 0:dx:a;
x2 = a+dx:dx:L;
C = 1;
    
dz1 = -F.*L.*b.^2./(6.*E.*I).*(3.*a./L.*(x1./L).^2-(1+2.*a./L).*(x1./L).^3);
dz2 = -F.*L.*a.^2./(6.*E.*I).*(3.*b./L.^3.*(L-x2).^2-(1+2.*b./L).*((L-x2).^3./L.^3));
D1 = ones(1,length(x1)).*F.*(b./L).^2.*(1+2.*a./L);
D2  = ones(1,length(x2)).*-F.*(a./L).^2.*(+2.*b./L);
MR1 = F.*b.^2.*(x1.*L+2.*x1.*a-a.*L)./(L.^3);
MR2 = F.*(b.^2.*x2.*L+2.*b.^2.*x2.*a-b.^2.*a.*L-L.^3.*x2+L.^3.*a)./(L.^3);
sigma1 = F.*u.*b.^2.*(x1.*L+2.*x1.*a-a.*L)./(I.*L.^3);
sigma2 = F.*u.*(b.^2.*x2.*L+2.*b.^2.*x2.*a-b.^2.*a.*L-L.^3.*x2+L.^3.*a)./(I.*L.^3);

dz = [dz1 dz2];
D = [D1 D2];
MR = [MR1 MR2];
sigma = abs([sigma1 sigma2]);

if a<b
    dz_max = -(2/3)*(F*a^2*b^3/(E*I*L^2))*(1/(1+(2*b/L)))^2;
    xdz = L-L*(1/(1+L/(2*b)));
    MR_max = (-F)*a*b^2/L^2;
    xMR = 0;
    sigma_max = abs((-F)*u*a*b^2/L^2/I);
    xsigma = 0;
elseif a>=b
    dz_max = -(2/3)*(F*a^3*b^2/(E*I*L^2))*(1/(1+(2*a/L)))^2;
    xdz = L*(1/(1+L/(2*a)));
    MR_max = (-F)*b*a^2/L^2;
    xMR = L;
    sigma_max = abs((-F)*u*b*a^2/L^2/I);
    xsigma = L;
end
Cz = abs(-3*E*I*L^3/(b^2*a^3*(L-a)));
xC = a;
end

%% Load case 11
if loadcase == 11;
C = 2;

dz = -q.*L.^4./(24.*E.*I).*((x./L).^2-2.*(x./L).^3+(x./L).^4);
D  = q.*L./2-q.*x;
MR = q.*(-6.*x.^2+6.*L.*x-L.^2)./12;
sigma  = abs(q.*u.*(-6.*x.^2+6.*L.*x-L.^2)./(12.*I));

dz_max = -q*L^4/(384*E*I);
xdz = L/2;
MR_max = -q*L^2/12;
xMR = [0 L];
sigma_max = abs(q*L^2*u/(12*I));
xsigma = [0 L];
Cz = abs(-384*E*I/L^4);
xC = L/2;
end

%% Load case 12
if loadcase == 12;
C = 1;

dz = -F.*x.^2.*(3.*L-2.*x)./(12.*E.*I);
D  = ones(1,length(x)).*F;
MR = (F.*L-2.*F.*x)./2;
sigma  = abs(u.*(F.*L-2.*F.*x)./(2.*I));

dz_max = -F*L^3/(12*E*I);
xdz = L;
MR_max = F*L/2;
xMR = [0 L];
sigma_max = abs(F*L*u/(2*I));
xsigma = [0 L];
Cz = abs(-12*E*I/L^3);
xC = L;

end

%% plots
scrsz = get(0,'ScreenSize');
pcs = 15;
formatSpec = '%6.4g';
figure('Position',[scrsz(1)+pcs scrsz(2)+4*pcs scrsz(3)-2*pcs scrsz(4)-10*pcs])
subplot 411; plot(x*1e3,dz*1e3,'r',[0,x(end)]*1e3,[0,0],'k:'); ylabel('\deltaz [mm]');
if C == 1
    title(['Loadcase ',num2str(loadcase),'       \deltax_m_a_x = ',num2str(dz_max*1e3),' [mm] @ x = ',num2str(xdz*1e3),' [mm],       |C_x_F| = ',num2str(Cz,formatSpec),' [N/m] @x_F = ',num2str(xC*1e3),' [mm]'])
elseif C == 2
    title(['Loadcase ',num2str(loadcase),'       \deltax_m_a_x = ',num2str(dz_max*1e3),' [mm] @ x = ',num2str(xdz*1e3),' [mm],       |C_x_\delta_x| = ',num2str(Cz,formatSpec),' [N/m^2] @x_\delta_x = ',num2str(xC*1e3),' [mm]'])
elseif C == 3
    title(['Loadcase ',num2str(loadcase),'       \deltax_m_a_x = ',num2str(dz_max*1e3),' [mm] @ x = ',num2str(xdz*1e3),' [mm],       |C_x_M| = ',num2str(Cz,formatSpec),' [Nm/rad] @x_M = ',num2str(xC*1e3),' [mm]'])
elseif C == 4
    title(['Loadcase ',num2str(loadcase),'       \deltax_m_a_x = ',num2str(dz_max*1e3),' [mm] @ x = ',num2str(xdz*1e3),' [mm],       |C_x_\delta_x| = ',num2str(Cz1,formatSpec),' [N/m^2] @x_\delta_x_1 = ',num2str(xC1*1e3),' [mm],      C_x_\delta_x = ',num2str(Cz2,formatSpec),' [N/m^2] @x_\delta_x_2 = ',num2str(xC2*1e3),' [mm]'])
end
subplot 412; plot(x*1e3,D,'r',[0,x(end)]*1e3,[0,0],'k:'); ylabel('F_s_h_e_a_r [N]');
subplot 413; plot(x*1e3,MR,'r',[0,x(end)]*1e3,[0,0],'k:'); ylabel('M_R_e_a_c_t_i_o_n [Nm]');
if xMR < 0
    title(['M_R_e_a_c_t_i_o_n_-_m_a_x = ',num2str(MR_max),' [Nm] @ x = Constant [mm]'])
elseif xMR >= 0
    title(['M_R_e_a_c_t_i_o_n_-_m_a_x = ',num2str(MR_max),' [Nm] @ x = ',num2str(xMR*1e3),' [mm]'])
end
subplot 414; plot(x*1e3,sigma/1e6,'r',[0,x(end)]*1e3,[0,0],'k:'); ylabel('|\sigma_b_e_n_d_i_n_g| [N/mm^2]'); xlabel('x [mm]: (0 - L)');
if xMR < 0
    title(['|\sigma_b_e_n_d_i_n_g_-_m_a_x| = ',num2str(sigma_max/1e6),' [N/mm^2] @ x = Constant [mm]'])
elseif xMR >= 0
    title(['|\sigma_b_e_n_d_i_n_g_-_m_a_x| = ',num2str(sigma_max/1e6),' [N/mm^2] @x = ',num2str(xsigma*1e3),' [mm]'])
end
