#!/usr/bin/octave
echo on
#Model bydunku pod trzesieniem ziemi 
#na i-te pietro dziala sila <--  k(i)(x(i) -x(i-1)) oraz k(i+1)(x(i+1)-x(i)) -->
#pierwsze pietro m_ix(1)''= - k(1)(x(2) -x(1)) 
#                m_ix(i)''= -k(i)(x(i) -x(i-1))+ k(i+1)(x(i+1)-x(i)) i=2... N-1
#ostatnie pietro m_ix(N)''= -k(N)(x(N) -x(N-1))
#
#W razie trzesienia na pierwsze pietro dziala 
#sila F(1)(t)=Gsin(a*t) a=3 typowe dla trzesienia
#

global A G a M #wspolczynniki 


function [A,B,M,K]=genmat(N,k=5000,m=10000,kk=k*ones(1,N),mm=m*ones(1,N))
#macierz 
 [m,n]=size(kk); if(m>n) kk=kk';endif
 M=diag(mm);
 K=diag(-kk)+diag(-[kk(2:N),0]) + diag(kk(2:N),1) + diag(kk(2:N),-1);
 A=M\K;
 B=[zeros(N,N), eye(N);A,zeros(N,N)];
endfunction

function y=pw(x,t=0)
#pole wektorowe 
#x''=AX+(G/mm(1))*sin(at)*E_1  E_1 piewszy wersor
#X1'=X2
#X2'=A*X1+ G*sin(at)*E_1
 global A G a M#wspolczynniki k m G a
 [m,N]=size(A);y=zeros(2*N,1);
 y(1:N)=x(N+1:2*N);
 y(N+1:2*N)=A*x(1:N);
 if(t<200)
   y(N+1)+=(G/M(1,1))*cos(a*t);
  else
    if(t<201)
      y(N+1)+=(G*(1-(t-200))/M(1,1))*cos(a*t);
    endif
 endif
endfunction


#sprawdzmy wartosci i wektory wlasne 
N=2;
k=5000#kg/s^2
m=10000; #kg
G=5000#kg/s^2 amplituda sily trzesienia 
a=3; #czestotliwosc sily wymuszajacej trzesienia ziemi

[A,B,M,K]=genmat(N);
lambda=eig(A);
#lambda(i)<0
#i*sqrt(-lambda(i)) - 
#omega=sqrt(-lambda(i))- i-ta czestotliwosc ukladu 
omega=sqrt(-lambda) #czestotliwosci budynku
T=2*pi./omega
#zalozmy ze uklad nie jest w stanie rownowagi tzn x_1(0)=randn
#ale pozostale x_k(0) i x_k'(0) = 0
#narysowac trajektorie dla t=[0,10] i [0,60]
G=0; #bez sily trzesienia ziemi
x0=zeros(2*N,1);
x0(1)=rand-0.5;
t=linspace(0,10);
X=lsode(@pw,x0,t);
Xs=X(:,1:N);max(abs(Xs)) #przesuniecia
#plot(t,Xs); #drgania wszystkich pieter
#drgania 1 i ostatniego
plot(t,Xs(:,[1,N]),[";1 i 2 p:k/m=.5;"]);
pause(2);
#dluzszy czas
t=linspace(0,30,1000);
X=lsode(@pw,x0,t);
Xs=X(:,1:N);max(abs(Xs)) #przesuniecia pieter
plot(t,Xs);pause(2); #drgania wszystkich pieter 

#to samo dla N=10
N=10;
k=5000#kg/s^2
m=10000; #kg
G=5000#kg/s^2 amplituda sily trzesienia 
a=3; #czestotliwosc sily wymuszajacej trzesienia ziemi
[A,B,M,K]=genmat(N);
lambda=eig(A);
#lambda(i)<0
#i*sqrt(-lambda(i)) - 
#omega=sqrt(-lambda(i))- i-ta czestotliwosc ukladu 
omega=sqrt(-lambda) #czestotliwosci budynku
T=2*pi./omega
#zalozmy ze uklad nie jest w stanie rownowagi tzn x_1(0)=randn
#ale pozostale x_k(0) i x_k'(0) = 0
#narysowac trajektorie dla t=[0,10] i [0,100]
G=0;
x0=zeros(2*N,1);
x0(1)=rand-0.5;
t=linspace(0,10);
X=lsode(@pw,x0,t);
Xs=X(:,1:N);max(abs(Xs)) #przesuniecia
#plot(t,Xs); #drgania wszystkich pieter dla N=10 za duzo
#drgania 1 i ostaniego
plot(t,Xs(:,[1,N]),[";1 i 10 p:k/m=.5;"]);
pause(2);
#dluzszy czas
t=linspace(0,30,1000);
X=lsode(@pw,x0,t);
Xs=X(:,1:N); #przesuniecia pieter
plot(t,Xs(:,[1,N]),[";1 i 10 p:k/m=.5;"]);pause(2); #drgania 2 i ostatniej kondygnacji




#trzesienie ziemi
G=5000;
x0=[zeros(2*N,1)];
t=linspace(0,10);
X=lsode(@pw,x0,t);
Xs=X(:,1:N);max(abs(Xs)) #przesuniecia
plot(t,Xs(:,[1,N])); pause(2);
t=linspace(0,60,1000);
X=lsode(@pw,x0,t);
Xs=X(:,1:N);Od=max(abs(Xs)) #przesuniecia
plot(t,Xs(:,[1,5,10]),[";1,5 i 10 p:k/m=.5;"]);pause(2); #drgania 1-szego 5-tego i ostatniego pietra


#zmienmy teraz a na czestotliwosc budynku 
a=2
X=lsode(@pw,x0,t);
Xs1=X(:,1:N);Od1=max(abs(Xs1))
plot(t,Xs(:,N),";trz o czest 3;",t,Xs1(:,N),";tr ziemi cz 2;");pause(2);
max(Od1./Od) #az razy wieksza amplituda....

#wezmy srednia z czestotliwosci wlasnych budynku
a=(omega(1)+omega(2))/2
X=lsode(@pw,x0,t);
Xs2=X(:,1:N);Od2=max(abs(Xs2))
plot(t,Xs(:,N),";trz o czest 3;",t,Xs2(:,N),";cz (om(1)+om(2)/)2;")
max(Od2./Od) 



#zwiekszmy G 2 krotnie i porwnajmy amplitude
G=10000; #G=2*G
X=lsode(@pw,x0,t);
Xs3=X(:,1:N);Od3=max(abs(Xs3))
plot(t,Xs(:,N),";trz o czest 3;",t,Xs3(:,N),";G podwojone;");
Od3./Od #czyli liniowo dla wszystkich pieter




#zmienmy budynek na k=10000 a m=5000 czyli 
N=10;
k=10000#kg/s^2
m=5000; #kg
G=5000#kg/s^2 amplituda sily trzesienia 
a=3; #czestotliwosc sily wymuszajacej trzesienia ziemi
[A,B,M,K]=genmat(N,k,m);
lambda=eigs(A);
#lambda(i)<0
omega=sqrt(-lambda) #czestotliwosci ukladu - tzn omega(i) -czestotliwosc drgania i-tego pietra
T=2*pi/omega
a=omega(1); #czestotliwosc trzesienia = czest naturalna budynku
t=linspace(0,10);
x0=[zeros(2*N,1)];
X=lsode(@pw,x0,t);
Xs=X(:,1:N); #przesuniecia
plot(t,Xs(:,[1,N]));pause(2);
t=linspace(0,60,1000);
X=lsode(@pw,x0,t);
Xs=X(:,1:N);Od=max(abs(Xs)) #max przesuniecia
plot(t,Xs(:,[1,N]));pause(2); #drgania 1 i ostatniego pietra pieter 


#a teraz znow standardowy ale z basenem na dachu tzn m_N=2*m
N=10;
k=5000#kg/s^2
m=10000; #kg
mm=m*ones(N,1);
kk=k*ones(N,1);
G=5000#kg/s^2 amplituda sily trzesienia 
a=3; #czestotliwosc sily wymuszajacej trzesienia ziemi
[A,B,M,K]=genmat(N,k,m,kk,mm);
lambda=eig(A);
#lambda(i)<0
#i*sqrt(-lambda(i)) - 
#omega=sqrt(-lambda(i))- i-ta czestotliwosc ukladu 
omega=sqrt(-lambda) #czestotliwosci budynku
T=2*pi./omega
x0=[zeros(2*N,1)];
t=linspace(0,40);
X=lsode(@pw,x0,t);
Xs=X(:,1:N);max(abs(Xs)) #przesuniecia
#a teraz z basenem
mm(1)*=2;
[A,B,M,K]=genmat(N,k,m,kk,mm);
lambda=eig(A);
omegab=sqrt(-lambda) #czestotliwosci budynku
Tb=2*pi./omegab
x0=[zeros(2*N,1)];
t=linspace(0,40);
X=lsode(@pw,x0,t);
Xsb=X(:,1:N);max(abs(Xsb)) 
plot(t,Xs(:,N),";std;",t,Xsb(:,N),";basen na 10p;");pause(2);
#odwrocona piramida
mm(N)=m;
for k=N-1:-1:1,
 mm(k)=mm(k+1)*.75;
endfor
#przeskalowujemy by laczna masa byla taka sama
mm=mm*(N*m)/sum(mm);
[A,B,M,K]=genmat(N,k,m,kk,mm);
lambda=eig(A);
omegab=sqrt(-lambda) #czestotliwosci budynku
Tb=2*pi./omegab
x0=[zeros(2*N,1)];
t=linspace(0,40);
X=lsode(@pw,x0,t);
Xso=X(:,1:N);max(abs(Xso)) 
for k=1:N,
 plot(t,Xs(:,k),";std;",t,Xso(:,k),[";",num2str(k),"pietro odw piram;"]);pause(2);
endfor


echo off
