clear all; close all

% Exercise 1 - a
t=linspace(0,100,1000);
k=16; m=1; c=0.1; wn=sqrt(k/m); xi=c/(2*sqrt(k*m));
wd=wn*sqrt(1-xi^2);
h=1/(m*wd)*exp(-xi*wn*t).*sin(wd*t);
figure;  plot(t,h); set(gca,'Fontsize',15);
xlabel('time'); ylabel('Impulse response h(t)')

% Exercise 1- b
w=linspace(0,20,10000);
H=1./(m*(wn^2-w.^2+2*j*xi*w*wn));
figure; set(gca,'Fontsize',15);
subplot(2,1,1); plot(w,20*log10(abs(H)));
xlabel('\omega');ylabel('|X/F|(dB)'); set(gca,'Fontsize',15);
subplot(2,1,2); plot(w,angle(H)*180/pi);
xlabel('\omega');ylabel('Phase(deg)'); set(gca,'Fontsize',15);

% Exercise 1- c
t=linspace(0,100,1000);
k=16; m=1; c=0.1; wn=sqrt(k/m); ci=[0.1 0.5 10];
color='rgb';

figure; 
for i=1:length(ci)
c=ci(i);
xi=c/(2*sqrt(k*m)); wd=wn*sqrt(1-xi^2);
h=1/(m*wd)*exp(-xi*wn*t).*sin(wd*t);
plot(t,h,color(i)); hold on;
end
set(gca,'Fontsize',15)
xlabel('time'); ylabel('Impulse response h(t)')

figure; 
for i=1:length(ci)
c=ci(i);
xi(i)=c/(2*sqrt(k*m)); wd=wn*sqrt(1-xi(i)^2);
H=1./(m*(wn^2-w.^2+2*j*xi(i)*w*wn));
subplot(2,1,1); plot(w,20*log10(abs(H)),color(i)); hold on
set(gca,'Fontsize',15)
xlabel('\omega');ylabel('|X/F|(dB)')
subplot(2,1,2); plot(w,angle(H)*180/pi,color(i)); hold on
xlabel('\omega');ylabel('Phase(deg)')
set(gca,'Fontsize',15)
end

% Values of xi
xi

%% Exercise 2

% Build impulse response
k=16; m=1; c=0.1;
xi=c/(2*sqrt(k*m));
wn=sqrt(k/m); wd=wn*sqrt(1-xi^2);
t=linspace(0,500*2*pi,20000);
h=exp(-xi*wn*t)/(m*wd).*sin(wd*t);

figure; plot(t,h); xlabel('t'); ylabel('h(t)'); set(gca,'Fontsize',15); set(gca,'XLim',[0 150])

% Sinusoidal inputs
% Ex 2-a
w=1*wn; f=sin(w*t); u=conv(h,f)*t(2); u=u(1:length(t)); 
figure; plot(t,f); xlabel('t'); ylabel('f(t)'); set(gca,'Fontsize',15); set(gca,'XLim',[0 150])
figure; plot(t,u); xlabel('t'); ylabel('x(t)'); set(gca,'Fontsize',15); set(gca,'XLim',[0 150])

% Ex 2-b
w=.95*wn; f=sin(w*t); u=conv(h,f)*t(2); u=u(1:length(t)); 
figure; plot(t,f); xlabel('t'); ylabel('f(t)'); set(gca,'Fontsize',15); set(gca,'XLim',[0 150])
figure; plot(t,u); xlabel('t'); ylabel('x(t)'); set(gca,'Fontsize',15); set(gca,'XLim',[0 150])

% Ex 2-c Random input
f=randn(length(t),1); 
u=conv(h,f)*t(2); u=u(1:length(t)); 
figure; plot(t,f); xlabel('t'); ylabel('f(t)'); set(gca,'Fontsize',15); set(gca,'XLim',[0 150])
figure; plot(t,u); xlabel('t'); ylabel('x(t)'); set(gca,'Fontsize',15); set(gca,'XLim',[0 150])

%%  Exercise 3

% Ex 3-a
load inputs
in(:,2:9)=in(:,2:9)*9.81; % to convert from g to m/s^2
figure; plot(in(:,1),in(:,2:end)); set(gca,'Fontsize',15) % plots all signals
xlabel('time(s)'); ylabel('acceleration(m/s^2)')
t=in(:,1); in=in(:,2:end);
frq=linspace(0,1/t(2),length(t));

% Ex 3-b Check frq content of inputs using FFT
figure;
for i=1:size(in,2)
E=t(2)*fft(in(:,i)); % need to multiply by dt to get FT.
plot(frq,abs(E)); hold on
end
set(gca,'XLim',[0 50]); set(gca,'Fontsize',15)
xlabel('Freq(Hz)'); ylabel('Amplitude of FFT')

% Ex 3-c build response spectra for each earthquake
% loop on k values
m=1; xi=0.05;
% Define period axis from 10-2 to 10, i.e. 100 Hz to 0.1 Hz
T0=logspace(-2,1,50);
% Determine k values to span the range of natural frequencies and
% corresponding wn values
k0=m./T0.^2*4*pi^2; w0=1./T0*2*pi; w0=w0';

% Compute response using Duhamel's integral and find max (Su value)
for i=1:length(k0)
k=k0(i); 
c=2*xi*sqrt(k*m); % for info only
wn=sqrt(k/m); wd=wn*sqrt(1-xi^2);
h=exp(-xi*wn*t)/(m*wd).*sin(wd*t);

for j=1:size(in,2)
 u=conv(h,in(:,j))*t(2); u=u(1:length(t));
 Su(i,j)=max(u);
end

end

figure; plot(T0,Su); set(gca,'Fontsize',15); xlabel('Period(s)'); ylabel('S_u(m)')

% Ex 3-d Pseudo acceleration spectrum
Sa=Su;

% multiply by wn^2
for i=1:size(Su,2)
    Sa(:,i)=w0.^2.*Su(:,i);
end

figure; plot(T0,Sa);
xlabel('Period(s)'); ylabel('Se(m/s^2)'); set(gca,'Fontsize',15)
set(gca,'XLim',[0 3.5])

% Mean spectrum
Sam=mean(Sa')'; figure; plot(T0,Sam); set(gca,'Fontsize',15)
xlabel('Period(s)'); ylabel('Se(m/s^2)')
set(gca,'XLim',[0 3.5])

% TB abd TC approximately 0.1s and 0.25s, i.e between 4 and 10 Hz


