%% Izračun notranjih napetosti in upogibnih momentov ravnega nosilca
%% Definicija količin
% Palica naj bo obremenjena tako s točkovno določenimi silami, točkovno podanimi navori in zvezno porazdeljenimi obremenitvami. Za numerično računanje moramo zvezne porazdelitve diskretizirati.

l=10; % dolžina v metrih
dx=0.1; % diskretizacija

x=0:dx:l; % vektor lege
n=length(x); % število diskretnih točk
%Za začetek postavimo vse na 0
fx=zeros(1,n); % točkovne sile
fy=zeros(1,n);
qx=zeros(1,n); % zvezno porazdeljene sile
qy=zeros(1,n);
m=zeros(1,n); % točkovni navori

% od polovice naprej naj bo palica obremenjena s 5kN/m
qy(round(n/2):n)=-5; % kaže navzdol
% na 7m od prijamališča naj bo diskretna sila 15kN
fy(round(7.0/10*n))=-15;

%% Izračun sil in navorov v podporah
% Sile in navore izračunamo tako, da uravnovesijo vse zunanje sile in navore, ki delujejo na palico. Tak sistem je v splošnem statično nedoločen, ker imamo šest neznank (sili v obeh smereh in navor, vsi v obeh krajiščih) in le tri enačbe (ravnovesje sil po obeh smereh in ravnovesje navorov). Sistem je statično določen le, če so nosilci zadosti "ohlapni" - vsaj tri od šestih navedenih količin morajo biti enake 0 (ali kakšni drugi vrednosti). Ne glede na vrsto opor, bodo rezultante zunanjih bremen brez opor enake.

%Najprej izračunamo zunanje sile brez opor
gx=0;
gy=0;
gm=0;
for j=1:n
    gx=gx+qx(j)*dx+fx(j);
    gy=gy+qy(j)*dx+fy(j);
    gm=gm+(qy(j)*dx+fy(j))*x(j)+m(j);
end

% V našem konkretnem primeru bomo upoštevali, kot da je palica na levem krajišču (A) nepomično vpeta, na desnem (B) pa pomično. Odtod sledi:
% F_{Ax}+F_x=0\\
% F_{Ay}+F_{By}+F_y=0\\
% F_{By}L+M=0

% Levo nepomična podpora, desno pomična
fax=-gx;
fb=-gm/l;
fay=-fb-gy;
% dodamo k točkovnim silam
fx(1)=fx(1)+fax;
fy(1)=fy(1)+fay;
fy(n)=fy(n)+fb;

%% Izračun notranjih napetosti in upogibnega momenta
% Ko enkrat imamo izračunane sile v podporah, lahko izračunamo še notranje napestosti in upogibni moment.

h=zeros(1,n);
v=zeros(1,n);
mn=zeros(1,n);
h(1)=-fx(1)-qx(1)*dx;
v(1)=-fy(1)-qy(1)*dx;
mn(1)=-v(1)*x(1)-m(1)-fy(1)*x(1)-qy(1)*x(1)*dx;
for j=2:n
    h(j)=h(j-1)-fx(j)-qx(j)*dx;
    t=fy(j)+qy(j)*dx; % Vmesni rezultat, ki ga bomo večkrat potrebovali
    v(j)=v(j-1)-t;
    mn(j)=mn(j-1)-(v(j)*x(j)-v(j-1)*x(j-1)+m(j)+t*x(j));
end
plot(x,h,'r',x,v,'b',x,mn,'g')

%