0==0;

#earthquake - model of n storejy building

#write a function computing the following 3diag matrix A=M\S
#for S 3diag 
# i-th row  k(i) -k(i)-k(i+1) k(i+1) i<n
#n-th               k(n) -k(n) i.e. k(n+1)=0
#M diagonal diag(m(1),...,m(n)]
function A=eqparam(k,m,n=length(m))
 k(n+1)=0;
 [kn,km]=size(k);
 if(kn>km) k=k'; end
 I=[1:n,2:n,1:n-1];
 J=[1:n,1:n-1,2:n];
 V=[-k(1:n)-k(2:n+1),k(2:n),k(2:n)];
 A=sparse(I,J,V,n,n);
 im=1.0./m;
 A=diag(im)*A;
end



function [omega,T]=freqper(A)
#computes frequencies and periods of building
 lambda=eig(A);
 omega=sqrt(-lambda);
 T=2*pi./omega;
end



function y=f(x,t)
global A m G gamma c 
 n=length(m);N=2*n;
 y=zeros(N,1);
 y(1:n)=x(n+1:N);
 y(n+1:N)=A*x(1:n) + c.*x(n+1:N);
 y(N)+=G*cos(gamma*t);
endfunction




#global parameters
global m k A G gamma c

gamma=3;
G=0; #no eathquake
n=2;
c=zeros(n,1); #no damping

k=m=ones(n,1);
k=5000*k;
m=10000*m;
A=eqparam(k,m); 
[omega,T]=freqper(A)
Tmax=10;N=101;h=Tmax/(N-1);
t=linspace(0,Tmax,N);
x0=dx0=zeros(n,1);x0(n)=1e-1;
X=lsode(@f,[x0;dx0],t);
figure(1);plot(t,X(:,1:2)',";no eq;");
amp=max(X(:,1:2));
mamp=max(amp);
xr=linspace(-5*mamp,5*mamp);
figure(2);
for j=1:2:length(t),
 td=["*g;t=" num2str(j*h) ";"];
 plot(xr,0*xr-0.1,"b",X(j,1:n),0:n-1,td,xr,0*xr+n,"b")
 pause(0.005);
end
figure(1);

#real earthquake
G=0.1
x0=dx0; #zero initial conditions
X=lsode(@f,[x0;dx0],t);
figure(1);plot(t,X(:,1:2)',";eaq n=2;");
amp=max(X(:,1:2));
mamp=max(amp);
xr=linspace(-5*mamp,5*mamp);
figure(2);
title("earthquake");
for j=1:2:length(t),
 td=["*g;t=" num2str(j*h) ";"];
 plot(xr,0*xr-0.1,"b",X(j,1:n),0:n-1,td,xr,0*xr+n,"b")
 pause(0.005);
end
figure(1);


#same 10 storey bldg
#global parameters
global m k A G gamma

gamma=3;
G=0; #no eathquake
n=10;
c=zeros(n,1); #no damping
k=m=ones(n,1);
k=5000*k;
m=10000*m;
A=eqparam(k,m); 
[omega,T]=freqper(A)
Tmax=10;N=101;h=Tmax/(N-1);
t=linspace(0,Tmax,N);
x0=dx0=zeros(n,1);x0(1)=1e-1;
X=lsode(@f,[x0;dx0],t);
plot(t,X(:,[1,n])',";no eaq n=10;");pause(1);
amp=max(X(:,1:n));
mamp=max(amp);
xr=linspace(-5*mamp,5*mamp);
figure(2);
for j=1:length(t),
 td=["*g;t=" num2str(j*h) ";"];
 plot(xr,0*xr-0.1,"b",X(j,1:n),0:n-1,td,xr,0*xr+n,"b")
 pause(0.005);
end
figure(1);

#to samo ale dla G=0.1; x0=dx0=0

#same 10 storey bldg
#global parameters
global m k A G gamma
gamma=3;
G=0.1; #earthquake
x0=dx0=zeros(n,1);
X=lsode(@f,[x0;dx0],t);
plot(t,X(:,[1,n])',";eaq n=10;");pause(1);
amp=max(X(:,1:n));
mamp=max(amp);
xr=linspace(-5*mamp,5*mamp);
figure(2);
for j=1:length(t),
 td=["*g;t=" num2str(j*h) ";"];
 plot(xr,0*xr-0.1,"b",X(j,1:n),0:n-1,td,xr,0*xr+n,"b")
 pause(0.005);
end
figure(1);

