%Vstupní hodnoty:
m=200; % [kg] Hmotnost motoru pro stojan (1/2 z celkové hmotnosti)
g=10; % [m*s^-2] Tíhové zrychlení
M=350; % [Nm] Zatížení od točivého momentu motoru pro stojan (1/2 z uvažovaného mometu)
F_nev=3000; % [N] Zatížení od nevyváženosti motoru pro stojan (1/2 z uvažované síly)
H=0.705; % [m] Výška stojanu
L=0.95; % [m] Šířka stojanu (vzdálenost mezi sloupy)
x=0; % [m] Působiště síly (0 -> střed sojanu)

%Vzniklé zatížení ze zadanách hodnot
F=(m*g)+F_nev;

%Matice
A=[(2/3)*H^3        0                                                                                               -H^2;
   0                ((1/3)*x^3)+((1/3)*((L/2)-x)^3)+(x*((L/2)-x)^2)+((x^2)*((L/2)-x))+((1/24)*L^3)+((1/2)*H*L^2)    -((1/2)*x^2)-((1/2)*((L/2)-x)^2)-(x*((L/2)-x))+((1/8)*L^2);
   -H^2             -((1/2)*x^2)-(x*((L/2)-x))-((1/2)*((L/2)-x)^2)+((1/8)*L^2)                                      x+((L/2)-x)+(2*H)+((1/2)*L)];
    
%Vektor neznámých 
NC=0;
TC=0;
MC=0;

X=[NC; TC; MC];

%Vektor výsledných sil
b=[-(1/2)*F*((L/2)-x)*H^2-(1/2)*M*H^2;
    -((1/3)*F*((L/2)-x)^3)-((1/2)*F*x*((L/2)-x)^2)-((1/4)*F*(L^2)*H)+((1/2)*F*L*x*H-((1/2)*M*((L/2)-x)^2)-(M*x*((L/2)-x))-((1/2)*M*L*H));
    ((1/2)*F*((L/2)-x)^2)+((1/2)*F*H*L)-(F*x*H)+(M*(L/2))-(M*x)+(M*H)];

%Výpočet
X=inv(A)*b;

NC=X(1,1);
TC=X(2,1);
MC=X(3,1);
%Reakce ve vetknutí
FAX=NC;
FAY=-TC-F;
MA=MC-(TC*(L/2))-(NC*H)-(F*((L/2)-x))-M;

FBX=-NC;
FBY=TC;
MB=MC+(TC*(L/2))-(NC*H);

%Moment v rohu nosníku
MD=MC-(F*(L/2))-(TC*((L/2)+x))-M;
ME=MC+(TC*(L/2));

%Výsledky vypsání
disp("Zadané parametry:");
fprintf("m="), disp(m);
fprintf("M="), disp(M);
fprintf("F_nev="), disp(F_nev);
fprintf("H="), disp(H);
fprintf("L="), disp(L);
fprintf("x="), disp(x);
disp("______________Reakce v místě C___________________");
fprintf("NC="), disp(NC);
fprintf("TC="), disp(TC);
fprintf("MC="), disp(MC);
disp("______________Reakce v místě A___________________");
fprintf("FAX="), disp(FAX);
fprintf("FAY="), disp(FAY);
fprintf("MA="), disp(MA);
disp("______________Reakce v místě B___________________");
fprintf("FBX="), disp(FBX);
fprintf("FBY="), disp(FBY);
fprintf("MB="), disp(MB);
disp("_____________Ohybový moment v rohu stojanu____________________");
fprintf("MD="), disp(MD);
fprintf("ME="), disp(ME);

%Hodnoty pro vykresleni grafů
MA_MAT=[];
MB_MAT=[];
MC_MAT=[];
MD_MAT=[];
ME_MAT=[];
FAy_MAT=[];
FAx_MAT=[];
FBy_MAT=[];
FBx_MAT=[];
poloha=[];
for x=0:0.005:((L/2)-0.225)
    poloha=[poloha x];
    A=[(2/3)*H^3        0                                                                                               -H^2;
         0                ((1/3)*x^3)+((1/3)*((L/2)-x)^3)+(x*((L/2)-x)^2)+((x^2)*((L/2)-x))+((1/24)*L^3)+((1/2)*H*L^2)    -((1/2)*x^2)-((1/2)*((L/2)-x)^2)-(x*((L/2)-x))+((1/8)*L^2);
      -H^2             -((1/2)*x^2)-(x*((L/2)-x))-((1/2)*((L/2)-x)^2)+((1/8)*L^2)                                      x+((L/2)-x)+(2*H)+((1/2)*L)];
    b=[-(1/2)*F*((L/2)-x)*H^2-(1/2)*M*H^2;
    -((1/3)*F*((L/2)-x)^3)-((1/2)*F*x*((L/2)-x)^2)-((1/4)*F*(L^2)*H)+((1/2)*F*L*x*H-((1/2)*M*((L/2)-x)^2)-(M*x*((L/2)-x))-((1/2)*M*L*H));
    ((1/2)*F*((L/2)-x)^2)+((1/2)*F*H*L)-(F*x*H)+(M*(L/2))-(M*x)+(M*H)];
    X=inv(A)*b;
    NC=X(1,1);
    TC=X(2,1);
    MC=X(3,1);
    MA_mat=MC-(TC*(L/2))-(NC*H)-(F*((L/2)-x))-M;
    MA_MAT=[MA_MAT MA_mat];
    MB_mat=MC+(TC*(L/2))-(NC*H);
    MB_MAT=[MB_MAT MB_mat];
    MC_mat=MC;
    MC_MAT=[MC_MAT MC_mat];
    MD_mat=MC-(F*((L/2)))-(TC*((L/2)+x))-M;
    MD_MAT=[MD_MAT MD_mat];
    ME_mat=MC+(TC*(L/2));
    ME_MAT=[ME_MAT ME_mat];
    FAy_mat=-TC-F;
    FAy_MAT=[FAy_MAT FAy_mat];
    FAx_mat=NC;
    FAx_MAT=[FAx_MAT FAx_mat];
    FBy_mat=TC;
    FBy_MAT=[FBy_MAT FBy_mat];
    FBx_mat=-NC;
    FBx_MAT=[FBx_MAT FBx_mat];
end

%Graf ohybových momentů
figure(1)
plot(poloha,MA_MAT,'LineWidth',3);
hold on;
plot(poloha,MB_MAT,'LineWidth',3);
plot(poloha,MC_MAT,'LineWidth',3);
plot(poloha,MD_MAT,'LineWidth',3);
plot(poloha,ME_MAT,'LineWidth',3);
hold off
title('Průběh ohybových momentů');
xlabel('Působiště zatížení x [m]');
ylabel('Ohybový moment [Nm]');
grid on
legend('MA','MB','MC','MD','ME');

%Graf reakčních sil
figure(2)
plot(poloha,FAy_MAT,'LineWidth',3);
hold on;
plot(poloha,FBy_MAT,'LineWidth',3);
plot(poloha,FAx_MAT,'LineWidth',3);
plot(poloha,FBx_MAT,'LineWidth',3);
hold off
title('Průběh reakčních sil');
xlabel('Působiště zatížení x [m]');
ylabel('Reakční síla [N]');
grid on
legend('FAy','FBy','FAx','FBx');

FAx_max=max(FAx_MAT); fprintf("FAx_max="), disp(FAx_max);
FAx_min=min(FAx_MAT); fprintf("FAx_min="), disp(FAx_min);

FAy_max=max(FAy_MAT); fprintf("FAy_max="), disp(FAy_max);
FAy_min=min(FAy_MAT); fprintf("FAy_min="), disp(FAy_min);

FBx_max=max(FBx_MAT); fprintf("FBx_max="), disp(FBx_max);
FBx_min=min(FBx_MAT); fprintf("FBx_min="), disp(FBx_min);

FBy_max=max(FBy_MAT); fprintf("FBy_max="), disp(FBy_max);
FBy_min=min(FBy_MAT); fprintf("FBy_min="), disp(FBy_min);

MA_max=max(MA_MAT); fprintf("MA_max="), disp(MA_max);
MA_min=min(MA_MAT); fprintf("MA_min="), disp(MA_min);

MB_max=max(MB_MAT); fprintf("MB_max="), disp(MB_max);
MB_min=min(MB_MAT); fprintf("MB_min="), disp(MB_min);

MC_max=max(MC_MAT); fprintf("MC_max="), disp(MC_max);
MC_min=min(MC_MAT); fprintf("MC_min="), disp(MC_min);

MD_max=max(MD_MAT); fprintf("MD_max="), disp(MD_max);
MD_min=min(MD_MAT); fprintf("MD_min="), disp(MD_min);

ME_max=max(ME_MAT); fprintf("ME_max="), disp(ME_max);
ME_min=min(ME_MAT); fprintf("ME_min="), disp(ME_min);
