%% ENGN0040 Homework 6 2020
%% Problem 1
%

clear all
m = 40; k = 38.6e03; c = 0.2518e03;
wn = sqrt(k/m);
z = c/(2*sqrt(k*m))
K = 1/k;
w = 100;
X0 = 1.e-03;
M = 1/sqrt((1-w^2/wn^2)^2 + (2*z*w/wn)^2);
F0 = X0/(K*M)
X0_resonance = K*F0/(2*z)

%
%% Problem 4
%

syms w z lam
M2 = ( (1-lam^2*w^2)^2 + (2*z*w)^2)/ ( (1-w^2)^2 + (2*z*w)^2);
wmax = solve(diff(M2,w)==0,w)
simplify(subs(M2,w,wmax))

%

clear all
syms w c z
m1 = 345.37;   m2 = 23.65;
L1 = 1; L2 = 9-L1;
k = 160.e06;
wn = sqrt(k*L1^2/(L1^2*m1 + m2*(L1+L2)^2));
lam = sqrt(m2*(L1+L2)*L2/(L1^2*m1 + (L1+L2)^2*m2));
z = c*L1/(2*sqrt(k*(m1*L1^2 + m2*(L1+L2)^2)));
M = sqrt( (1-lam^2*w^2/wn^2)^2 + (2*z*w/wn)^2)/ sqrt( (1-w^2/wn^2)^2 + (2*z*w/wn)^2);
watmax = eval(solve(diff(M,w)==0,w));
Matmax = simplify(subs(M,w,watmax(2)));
csol = vpasolve(Matmax==2,c,[1000,10e06])
figure
fplot(subs(M,c,csol),[0,4*wn],'Color',[1 0 0],'LineWidth',2,'DisplayName','\zeta=0.1')
ylabel('Magnification');
xlabel({'Frequency (rad/s)'});
axes1 = gca;
ylim([0 2.5])
box(axes1,'on');
set(axes1,'FontName','Times New Roman','FontSize',16,'XGrid','on','YGrid',...
    'on');
eval(subs(M,[w,c],[wn,csol]))