clear
clc
close all
M = 3.98*10^7;
m = 6.6*10^5;
k_f = 35*10^6;
k_d1 = 580.4*10^3;
k_d2 = 546.1*10^3;
c_d = 37.135*10^3;
F_0 = 50*10^3;
omega = 0.7:0.001:1.2;
omega_1 = sqrt(k_f/M);
omega_2_1 = sqrt(k_d1/m);
omega_2_2 = sqrt(k_d2/m);
r = omega/omega_1;
alpha1 = omega_2_1/omega_1;
alpha2 = omega_2_2/omega_1;
mu = m/M;
delta_st = F_0/k_f;
zeta = c_d/(2*m*omega_1);
for i = 1:length(omega)
X_M1(i) = F_0/(sqrt((k_f-(M*omega(i)^2))^2));
end
X_m1 = 0;
for i = 1:length(omega)
X_M2(i) = abs((alpha1^2-r(i)^2)/(((1-r(i)^2)*(alpha1^2-r(i)^2))-(mu*alpha1^2*r(i)^2)))*delta_st;
end
for i = 1:length(omega)
X_m2(i) = abs((alpha1^2)/((1-r(i)^2)*(alpha1^2-r(i)^2)-(mu*alpha1^2*r(i)^2)))*delta_st;
end
for i = 1:length(omega)
X_M3(i) = sqrt(((alpha2^2-r(i)^2)^2+(2*zeta*r(i))^2)/((((1-r(i)^2)*(alpha2^2-r(i)^2))-(mu*alpha2^2*r(i)^2))^2+((2*zeta*r(i))*(1-r(i)^2-mu*r(i)^2))^2))*delta_st;
end
for i = 1:length(omega)
X_m3(i) = sqrt(((alpha2^2)^2+(2*zeta*r(i))^2)/(((1-r(i)^2)*(alpha2^2-r(i)^2)-(mu*alpha2^2*r(i)^2))^2+((2*zeta*r(i))*(1-r(i)^2-mu*r(i)^2))^2))*delta_st;
end
zeta_optimal = sqrt((mu*(3-(sqrt(mu/(mu+2)))))/(8*(1+mu)^3));
tuned_alpha = 1/(1+mu);
for i = 1:length(omega)
X_M4(i) = sqrt(((tuned_alpha^2-r(i)^2)^2+(2*zeta_optimal*r(i))^2)/((((1-r(i)^2)*(tuned_alpha^2-r(i)^2))-(mu*tuned_alpha^2*r(i)^2))^2+((2*zeta_optimal*r(i))*(1-r(i)^2-mu*r(i)^2))^2))*delta_st;
end
for i = 1:length(omega)
X_m4(i) = sqrt(((tuned_alpha^2)^2+(2*zeta_optimal*r(i))^2)/(((1-r(i)^2)*(tuned_alpha^2-r(i)^2)-(mu*tuned_alpha^2*r(i)^2))^2+((2*zeta_optimal*r(i))*(1-r(i)^2-mu*r(i)^2))^2))*delta_st;
end
subplot(2,1,1)
plot(r,X_M1,'LineWidth',2)
hold on
plot(r,X_M2,'-.r','LineWidth',2)
hold on
plot(r,X_M3,'.g','LineWidth',1)
hold on
plot(r,X_M4,'--b','LineWidth',1)
axis([0.5 1.5 0 0.05])
ylabel('X_M')
xlabel('r')
legend('Skyscraper w/ no mass damper','Skyscraper w/ undamped mass','Skyscraper w/ damped mass','Optimally Tuned')
subplot(2,1,2)
plot(r,X_m1,'LineWidth',2)
hold on
plot(r,X_m2,'-.r','LineWidth',2)
hold on
plot(r,X_m3,'.g','LineWidth',2)
hold on
plot(r,X_m4,'--b','LineWidth',1)
axis([0.5 1.5 0 0.25])
ylabel('X_m')
xlabel('r')
[pksM1,locsM1] = findpeaks(X_M1);
[pksM2,locsM2] = findpeaks(X_M2);
[pksM3,locsM3] = findpeaks(X_M3);
[pksM4,locsM4] = findpeaks(X_M4);