clear all;
close all;
n=5; %pocet podpernych tyci
m=19;
     %pocet lamel
a=150; %vzdalenost dvou tyci od sebe [mm]
c=150; %vzdalenost krajnich tyci od konce lamely [mm]
b=19.5; %vzdalenost os dvou sousednich lamel [mm]
e=3; %poloha zatezujici sily Q (pro pocet tyci 4, zatizeni v polovine odpovida e=2.5 
p=10; %cislo lamely, ktera je zatizena silou Q (pocitano od horniho okraje)
Q=981; %zatezujici osamela sila [N]
EL=72000; %modul pruznosti v tahu pro lamelu [MPa]
ET=210000; %modul pruznosti v tahu pro tyc [MPa]
JL=1742; %kvadriticky moment pruezu lamely k ose z [mm^4]
WoL=205;
WoT=50.2;
vyskaL=17;
JT=201; %kvadriticky moment prurezu lamely k ose z [mm^4], kruh 8mm->200 mm^4, ctverec 8x8-> 341 mm^4, 10x10->833 mm^4
vyskaT=10;
vel=2*m+2*n+m*n;
% vypocet konstant %
for v=1:1:n
    mL=(a*(n-v)+c)/((n-1)*a+2*c);
    r=((v-1)*a+c)/((n-1)*a+2*c);
    vc=(v-1)*a+c;
    nv=(n-v)*a+c;
    for u=1:1:n
        uc=(u-1)*a+c;
        if u<v
            ku=-mL*((vc^3-uc^3)/3-uc*(vc^2-uc^2)/2);
            kU(u,v)=ku/(EL*JL);
        end
        nu=(n-u)*a+c;
        if u>v
           ku=-r*((nv^3-nu^3)/3-nu*(nv^2-nu^2)/2);
           kU(u,v)=ku/(EL*JL);
       end
   end
   kl=-mL*(vc^3/3);
   kL(v)=kl/(EL*JL);
   kr=-r*(nv^3/3);
   kR(v)=kr/(EL*JL);
   ec=(e-1)*a+c;
   ne=(n-e)*a+c;
   if e>v
       kf1=-r*((nv^3-ne^3)/3-ne*(nv^2-ne^2)/2);
       kF(v)=kf1/(EL*JL);
   end
   if e<v
       kf2=-mL*((vc^3-ec^3)/3-ec*(vc^2-ec^2)/2);
       kF(v)=kf2/(EL*JL);
   end
end

for s=1:1:m
    h=(m+1-s)/(m+1);
    d=s/(m+1);
    ms=(m+1-s)*b;
    for t=1:1:m
        mt=(m+1-t)*b;
        if s>t
            kt=-h*((s^3-t^3)*b^3/3-t*b^3*(s^2-t^2)/2);
            kT(t,s)=kt/(ET*JT);
        end
        if s<t
            kt=-d*((ms^3-mt^3)/3-mt*(ms^2-mt^2)/2);
            kT(t,s)=kt/(ET*JT);
        end
    end
    kh=-h*(s*b)^3/3;
    kH(s)=kh/(ET*JT);
    kd=-d*(ms^3)/3;
    kD(s)=kd/(ET*JT);
end
% matice promennych 
A=zeros(vel);
% zaplneni matice
% silova rovnovaha, sily L, R, F(t, u)
for i=1:1:m
    jL=i;
    A(i,jL)=1;
    jR=i+m;
    A(i,jR)=1;
    for jF=2*m+n*(i+1)+1:1:2*m+n*(i+2);
        A(i,jF)=1;
    end
end
% momentove rovnice sil R,F(t,u)
for i=m+1:1:2*m
    j=i;
    A(i,j)=(n-1)*a+2*c;
end
for i=m+1:1:2*m
    k1=i-m;
    for j=2*m+n*(k1+1)+1:1:2*m+n*(k1+2);
        k2=j-(2*m+n*(k1+1));
        A(i,j)=(k2-1)*a+c;
    end
end
% silova rovnovaha sil H, D, F(t,u)
for i=2*m+1:1:2*m+n
    jH=i;
    A(i,jH)=1;
    jD=i+n;
    A(i,jD)=1;
end
for k=1:1:m
    for i=2*m+1:1:2*m+n
        j=i+n*(k+1);
        A(i,j)=1;
    end
    for im=2*m+n+1:1:2*m+2*n
        jm=im+n*k;
        A(im,jm)=k*b;
    end
end
%momentovka D
for i=2*m+n+1:1:2*m+2*n
    j=i-n;
    A(i,j)=(m+1)*b;
end
% def.podminky sily L
for j=1:1:m
    for i=2*m+n*(1+j)+1:1:2*m+n*(2+j)
        v=i-(2*m+n*(1+j));
        A(i,j)=kL(v);
    end
end
%def. podminky od sily R
for j=m+1:1:2*m
    for i=2*m+n*(1+j-m)+1:1:2*m+n*(2+j-m)
        v=i-(2*m+n*(1+j-m));
        A(i,j)=kR(v);
    end
end
%def. podminky sily H
for k=1:1:m
    for i=2*m+n*(k+1)+1:1:2*m+n*(k+2)
        j=i-n*(k+1);
        s=k;
        A(i,j)=kH(s);
    end
end
%def. podmnky D
for k=1:1:m
    for i=2*m+n*(k+1)+1:1:2*m+n*(k+2)
        j=i+n-n*(k+1);
        s=k;
        A(i,j)=kD(s);
    end
end
% def podminky F(t,u) od lamel
for t=1:1:m
    for v=1:1:n
        for u=1:1:n
            i=2*m+n*(1+t)+v;
            j=2*m+n*(1+t)+u;
            A(i,j)=A(i,j)+kU(u,v);
        end
    end
end
%def. podminky F(t,u) od tyci
for s=1:1:m
    for u=1:1:n
        for t=1:1:m
            i=2*m+n*(1+s)+u;
            j=2*m+n*(1+t)+u;
            A(i,j)=A(i,j)+kT(t,s);
        end
    end
end
A;
% tvorba vektoru b, prave strany
for i=1:1:2*m+2*n+(p-1)*n
    bb(i)=0;
    if i==p
        bb(i)=-Q;
    else
    end
    if i==m+p
        bb(i)=-Q*((e-1)*a+c);
    else
    end
end
for i=2*m+2*n+(p-1)*n+1:1:2*m+2*n+p*n
    v=i-(2*m+2*n+(p-1)*n);
    bb(i)=-Q*kF(v);
end
for i=2*m+2*n+p*n+1:1:vel
    bb(i)=0;
end
bv=bb'; %vektor prave strany
x=A\bv; %vektor reakcnich sil
%vektory sil
for i=1:1:m
    L(i)=x(i);
end
for i=m+1:1:2*m
    j=i-m;
    R(j)=x(i);
end
for i=2*m+1:1:2*m+n
    j=i-2*m;
    H(j)=x(i);
end
for i=2*m+n+1:1:2*m+2*n
    j=i-(2*m+n);
    D(j)=x(i);
end
for k=1:1:m
    for i=2*m+n*(1+k)+1:1:2*m+n*(2+k)
        j=i-(2*m+n*(1+k));
        F(k,j)=x(i);
    end
end
L;
R;
H;
D;
F;
% vypocet maximalnich ohybovych momentu
nejk=0;
for k=1:1:n
    nejk=k;
    if k>=e
        break
    end
end
nejk;
%cyklus pro lamelu
for i=1:1:m
    vektQ(i)=0;
end
vektQ(p)=Q;
vektQ;
%prvni cast
for t=1:1:m
    for i=1:1:nejk-1
        Moa=0;
        Mola=-L(t)*((i-1)*a+c);
        for u=1:1:i-1
            Moa=Moa-F(t,u)*(i-u)*a;
        end
        Mo(t,i)=Moa+Mola;
    end
%druha cast
    Mob=0;  
    for u=1:1:nejk-1
        Mob=Mob-F(t,u)*(e-u)*a;
    end
    Mo(t,nejk)=Mob-L(t)*((e-1)*a+c);
%treti cast
    Moc=0;
    for u=1:1:nejk-1
        Moc=Moc-F(t,u)*(nejk-u)*a;
    end
    Mo(t,nejk+1)=Moc-L(t)*((nejk-1)*a+c)-vektQ(t)*(nejk-e)*a;
%ctvrta cast
    for i=nejk+1:1:n
        Mod=0;
        for u=1:1:i-1
            Mod=Mod-F(t,u)*(i-u)*a;
        end
        Mo(t,i+1)=Mod-L(t)*((i-1)*a+c)-vektQ(t)*(i-e)*a;
    end
end
Mo;
maxL=0;
for i=1:1:m
    for j=1:1:n
        if abs(Mo(i,j))>maxL
            maxL=abs(Mo(i,j));
        end
        
    end
end
maxL;
%ohybovy moment v tycich
for u=1:1:n
    for i=1:1:m
        Mox=0;
        for t=1:1:i-1
            Mox=Mox+F(t,u)*(i-t)*b;
        end
        Motyce(u,i)=Mox+D(u)*i*b;
    end
end
Motyce;
maxT=0;
for i=1:1:n
    for j=1:1:m
        if abs(Motyce(i,j))>maxT
            maxT=abs(Motyce(i,j));
        end
        
    end
end
maxT;
%maximalni napeti v lamela
sigmaL=maxL/WoL
%maximalni napeti v tyci
sigmaT=maxT/WoT