% Dynamika Budowli - cwiczenie nr 5 (szablon do wyznaczania K i M)
% Autor: Magdalena Rucka
% rok. akad. 2014/15

clear;clc;close all

Ne  = 21          % liczba elementow
Neq =   37        % liczba niewiadomych uogolnionych przemieszczen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Dane materialowe i geometryczne
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x1 = 3;
x2 = 2.4;
x3 = 2.1;

E=  (20+x1)*10^9       ;    % modul sprezystosci [N/m2]
H1 = 4 + x2;
H2 = 4.5 + x2 ;   % wysokosc ramy [m]
L= 3 + x3;             % rozpietosc ramy [m]
b1 = (20 + x1)*10^-2  ;          % szerokosc przekroju belki [m]
h1=   (30 + x1)*10^-2  ;        % wysokosc przekroju belki [m]
b2=   (30 + x3)*10^-2   ;       % szerokosc przekroju slupa [m]
h2=   (35 + (x2 + x3))*10^-2  ;        % wysokosc przekroju slupa [m]
ro=   2500    ;      % gestosc [kg/m3]

Ib= (b1*h1^3)/12;            % moment bezwladnosci belki [m4]
Ic=  (b2*h2^3)/12;           % moment bezwladnosci slupa [m4]
mi_b=  ro*b1*h1;         % masa rozlozona belki [kg/m]
mi_c=  ro*b2*h2;         % masa rozlozona slupa [kg/m]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        TABLICA CECH ELEMENTOW
%-------------------------------------------
%        Ne     EI     L    ro
%-------------------------------------------
CE   = [ 1    E*Ic   H1/3   mi_c
         2    E*Ic   H1/3   mi_c
         3    E*Ic   H1/3   mi_c
         4    E*Ic   H2/3   mi_c
         5    E*Ic   H2/3   mi_c
         6    E*Ic   H2/3   mi_c
         7    E*Ib   L/3   mi_b
         8    E*Ib   L/3   mi_b
         9    E*Ib   L/3   mi_b
         10    E*Ic   H2/3   mi_c
         11    E*Ic   H2/3   mi_c
         12    E*Ic   H2/3   mi_c
         13    E*Ib   L/3   mi_b
         14    E*Ib   L/3   mi_b
         15    E*Ib   L/3   mi_b
         16    E*Ib   L/3   mi_b
         17    E*Ib   L/3   mi_b
         18    E*Ib   L/3   mi_b
         19    E*Ic   H1/3   mi_c
         20    E*Ic   H1/3   mi_c
         21    E*Ic   H1/3   mi_c];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        TABLICA WEKTOROW ALOKACJI
%-------------------------------------------
%        Ne | Vi | Fi | Vk | Fk
%-------------------------------------------
ALOK=[   1    0    1    2    3
         2    2    3    4    5
         3    4    5    6    7
         4    6    7    8    9
         5    8    9   10   11
         6   10   11   17   12
         7   0    12   13   14
         8   13   14   15   16
         9   15   16   0    18
         10  19   20   17   18
         11  21   22   19   20
         12  6    23   21   22
         13  0    7    24   25
         14  24   25   26   27
         15  26   27   0    23
         16  0    23   28   29
         17  28   29   30   31
         18  30   31   32   33
         19  34   35   6    23
         20  36   37   34   35
         21  0    0    36   37];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Agregacja globalnej macierzy sztywnosci i mas
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[K]=agreg_k_beam(Neq,Ne,CE,ALOK);
[M]=agreg_m_beam(Neq,Ne,CE,ALOK);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Rozwiazanie problemu wlasnego
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[mode,vale] = eig(K,M);
omega=sqrt(diag(vale));
f=omega/2/pi
[mode,f]=sort_norm(mode,f)
postac_1_2_3=mode(:,1:3)
 Sk=trace(K)
 Sm=trace(M)