% eigenproblem of 3-dof plane buildnig system
% with equal stiffness
%k=1.300*10^6 N/m & equal masses m=1000kg
clear;
clc;
%
ones=[1;1;1];
% zaczynamy od jednostkowego wektora
BB=ones'.*1000;
% masy poszczególnych piêter to zwyk³y wektor (po 1000kg)
B=diag(BB);
% ogólnie macierz bezw³adnosci jest diagonalna
% Deklaracja macierzzy sztywnoœci wierszami [K]*1.296*10^6
K=[2 -1 0; -1 2 -1; 0 -1 1].*1.296E6;
% sztywnoœæ
k=1.296*10^6
% obliczenie macierzy A do zagadnienia wlasnego w postaci standardowej A=lambda*q
A=B^(-1)*K;
% sprowadzenie do postaci standardowej problemu w³asnego
[V,D]=eig(A);
% rozwi¹zanie zagadnienia w³asnego macierzy A
% umieszczenie omega^2=lambda na diagonali w D oraz
% kolejnych wektorów w³asnych w macierz V
% umieszczenie kolejnych okresow drgan w wektorze T:
T(1)=2.*pi./sqrt(D(1,1));T(2)=2.*pi./sqrt(D(2,2));T(3)=2.*pi./sqrt(D(3,3));
T
% normalizacja wektorów w³asnych wzgl max, co do bezwz. wartoœci elementu
s=1.0;
for j=1:3
    max=0.0;
    for i=1:3
        if abs(V(i,j))>max
            max=abs(V(i,j));
            s=sign(V(i,j));
        end %if
    end %for i
    for i=1:3
        W(i,j)=s.*V(i,j)./max;
    end %for i
end %for j
% koniec normalizacji -> wyniki w w macierzy W
W(:,1) % wydruk kolumny 1 zaczynaj¹c od do³u konstrukcji
W(:,2) % wydruk kolumny 2 zaczynaj¹c od do³u konstrukcji
W(:,3) % wydruk kolumny 3 zaczynaj¹c od do³u konstrukcji
% wykresy form w³asnych
% figure;
subplot(2,2,1)
plot([0,W(1,3),W(2,3),W(3,3)],[0,1,2,3]);grid on;
title('Trzecia postaæ drgañ')
xlabel('Przemieszczenia')
ylabel('Wysokoœæ')
subplot(2,2,2)
plot([0,W(1,2),W(2,2),W(3,2)],[0,1,2,3]);grid on; 
title('Druga postaæ drgañ')
xlabel('Przemieszczenia')
ylabel('Wysokoœæ')
subplot(2,2,3)
plot([0,W(1,1),W(2,1),W(3,1)],[0,1,2,3]);grid on;
title('Pierwsza postaæ drgañ')
xlabel('Przemieszczenia')
ylabel('Wysokoœæ')
% sprawdzanie ortogonalnoœci otrzymanych wektorów w³asnych
% 'orthogonality with respect to B matrix'
% W(:,3)'*B*W(:,1) % orthogonality with inertia holds
% W(:,2)'*B*W(:,1) % orthogonality with inertia holds
% 'orthogonality with respect to B matrix'
% W(:,3)'*K*W(:,1) % orthogonality with inertia holds
% W(:,2)'*K*W(:,1) % orthogonality with inertia holds