%%
clear all; close all

%% Exercise 1

%% Ex1 a) Build matrices of the system
m=1; k=16; k1=16; c=0.1; c1=0.1;
K=[k1+k -k 0 ; -k 2*k -k ; 0 -k k];
M=[m 0 0 ; 0 m 0 ; 0 0 m];
C=[c1+c -c 0 ; -c 2*c -c ; 0 -c c];

%% Ex 1 b) compute and represent modeshapes
[V,D]=eig(K,M);
f=1/(2*pi)*sqrt(diag(D)); % frequencies in Hz
figure;set(gcf,'PaperUnits','centimeters','PaperPosition',[0 0 15 15],...
    'PaperSize',[15 15])
subplot(3,1,1); barh(V(:,1)); set(gca,'fontsize',15);title(['mode 1 / ' num2str(f(1)) ' Hz']);set(gca,'fontsize',15);
subplot(3,1,2); barh(V(:,2)); set(gca,'fontsize',15);title(['mode 2 / ' num2str(f(2)) ' Hz']);set(gca,'fontsize',15);
subplot(3,1,3); barh(V(:,3)); set(gca,'fontsize',15);title(['mode 3 / ' num2str(f(3)) ' Hz']);set(gca,'fontsize',15);

% Check orthogonality
Mr=V'*M*V % mass normalized
Kr=V'*K*V
D  % The stiffness matrix projected in modal basis is D


%% Ex 1 c- Projection in the modal basis
mui=diag(Mr); wi=sqrt(diag(Kr))./mui;
xi=(1./(2*mui.*wi)).*diag(V'*C*V); % modal damping
wd=wi.*(sqrt(1-xi.^2)); 
F=V'*[0;0;1]; % modal forces

% Impulse response for dof 3
t=linspace(0,100,1000);
u=zeros(1,length(t));
for i=1:3
    u=u+V(3,i)*(F(i)*exp(-xi(i)*wi(i)*t)/(mui(i)*wd(i))).*sin(wd(i)*t);
end
figure;
plot(t,u,'linewidth',2,'color','black'); set(gca,'fontsize',15);
xlabel('time (sec)'); ylabel('Amplitude (m)')
title('impulse response - x_3');


%% Ex 1 d) Bode diagram
w=linspace(0,10,1000);
U=zeros(1,length(w));

for i=1:3
 U=U+ (F(i)*V(3,i))./(w.^2 - wi(i)^2 + 2*j*xi(i)*wi(i)*w);
end

figure;
semilogy(w/(2*pi),abs(U),'k','linewidth',2); hold on;;set(gca,'fontsize',15)
xlabel('Freq(Hz)'); ylabel('Ampl X3/F (dB)')

xi(2)=xi(2)*5;
U2=zeros(1,length(w));
for i=1:3
 U2=U2+ (F(i)*V(3,i))./(w.^2 - wi(i)^2 + 2*j*xi(i)*wi(i)*w);
end
semilogy(w/(2*pi),abs(U2),'r','linewidth',2)
legend('initial','xi(2)*5')
title('Bode diagram - displacement x_3')
% increasing the damping of mode 2 leads to a reduction of the second peak only

% Acceleration
A=(-w.^2).*U;
semilogy(w/(2*pi),abs(A),'g','linewidth',2);set(gca,'fontsize',15)
xlabel('Freq(Hz)'); ylabel('Ampl (dB)')
title('Bode diagram - acceleration x_3')
legend('initial','xi(2)*5','initial acc','location','southeast')

%% Ex 1 e) c1=0;
c1=0;
C=[c1+c -c 0 ; -c 2*c -c ; 0 -c c];

% Direct reference method
U=zeros(3,length(w));
for i=1:length(w)
    U(:,i)=(K-w(i)^2*M + 1j*w(i)*C)\[0 0 1]';
end

figure;
semilogy(w/(2*pi),abs(U(3,:)),'k','linewidth',2); hold on ;set(gca,'fontsize',15)
xlabel('Freq(Hz)'); ylabel('Ampl (dB)')

% Modal damping hypothesis
xi=(1./(2*mui.*wi)).*diag(V'*C*V);
U2=zeros(1,length(w));
for i=1:3
 U2=U2+ (F(i)*V(3,i))./(w.^2 - wi(i)^2 + 2*1j*xi(i)*wi(i)*w);
end
semilogy(w/(2*pi),abs(U2),'r','linewidth',2)
legend('direct method','modal approach')
title('Bode diagram - displacement x_3')

% The modal damping hypthesis is valid, there is only a very minor
% differente at anti-resonances.

print(gcf,'-dpdf','BODE_B_U')

%% Ex 2

%% Ex 2 a) : Mass and stiffness matrices
% number of floors
n=5;

% Stiffness of floors :
k=1e9;
% mass of floors
m=7e4;

% Build stiffness and mass matrix (loop to assemble)
M=diag(m*ones(n,1));
K=zeros(n,n); K(1,1)=k;
for i=2:n
    K(i-1:i,i-1:i) = K(i-1:i,i-1:i) + [k -k; -k k];
end

%% Ex 2 b) Compute and represent modeshapes
[md,D]=eig(K,M);
f=sqrt(diag(D))/(2*pi); wn=2*pi*f;

figure; 
for i=1:n
subplot(1,n,i); barh(md(:,i)); xlabel(['mode ' num2str(i)]); set(gca,'Fontsize',15)
end

%% Ex 2 c) Compute response from 0 to 100 Hz
w=linspace(0,100,1000)*2*pi;

% Modal approach (mui=1 for all modes)
xi=0.01*ones(n,1);
U=zeros(1,length(w));
Gamma=-md'*M*ones(n,1); % unit acceleration of basis
for i=1:n
    U=U+ (Gamma(i))./(wn(i)^2-w.^2+2*1i*xi(i)*w*wn(i))*md(5,i);
end

figure; semilogy(w/(2*pi),abs(U),'linewidth',2);
xlabel('Freq(Hz)'); ylabel('$|\frac{X_{5r}}{-\omega^2 X_0}|$','Interpreter','latex'); set(gca,'Fontsize',15); set(gca,'XLim',[0 50])
a=get(gcf,'CurrentAxes'); set(a,'Box','off'); set(a,'Linewidth',2)

%% Ex 2 d) compute impulse response
mui=diag(md'*M*md); wi=wn;
wd=wi.*(sqrt(1-xi.^2)); 
F=Gamma; % modal forces

% Impulse response for dof 5 (for a unit ground acceleration in m/s^2)
t=linspace(0,20,1000);
u=zeros(1,length(t));
for i=1:3
    u=u+md(5,i)*(F(i)*exp(-xi(i)*wi(i)*t)/(mui(i)*wd(i))).*sin(wd(i)*t);
end
figure;
plot(t,u,'linewidth',2,'color','black'); set(gca,'fontsize',15);
xlabel('time (sec)'); ylabel('Amplitude (m)')

% load earthquake and make convolution
load inputs
f=in(:,2)*9.81; t=in(:,1); % it is in g , convert to m/s^2
u5r=t(2)*conv(f,u); u5r=u5r(1:length(t));

figure;
plot(t,u5r,'linewidth',2,'color','black'); set(gca,'fontsize',15);
xlabel('time (sec)'); ylabel('Amplitude (m)')
title('Relative displ fifth floor')

% use diff(:,2) to compute second derivative (length of res. vector is shorter)
a5r=diff(u5r,2); a5=a5r+f(1:length(a5r));

figure;
plot(t(1:length(a5r)),a5,'linewidth',2,'color','black'); set(gca,'fontsize',15);
xlabel('time (sec)'); ylabel('Amplitude (m/s^2)')
title('Absolute acceleration fifth floor')

