% AbsorberDesign.m   
% by K.C. Park  10  April 2005
clear all;

% Define some useful numbers:
dtr=pi/180;
Hz2rps=2*pi; rps2Hz=1/2/pi; % Conversions to/from Hz and rad/s

% Natural frequencies of individual masses: 10Hz
Omega1=10*Hz2rps;
mu = 0.01;
%
omega2=(1./(1+mu))*Omega1; % optimized frequency ratio
%
%  omega2 = 9.5*Hz2rps; 

zeta_opt = sqrt(3*mu/(8*(1+mu)^3));  

zeta_data = [0.5  1.0  2.0]*zeta_opt;

for i =1:3,  %zeta is the damping ratio
%
zeta = zeta_data(1,i);
% Define a system in transfer function form
a0=1; a1 =2*zeta*omega2;  a2 = omega2^2;
%
b0 = 1; 
b1 = 2*zeta*(1+mu)*Omega1; 
b2 = Omega1^2+ (1+mu)* omega2^2;
b3 = 2*zeta*Omega1^3; 
b4 = Omega1^2*omega2^2;

% develop transfer function
sys = tf ( [a0 a1 a2], [b0 b1  b2  b3 b4]);

% Pick a range of frequencies to look at from 0 Hz to 100 Hz
wf= logspace(0.5, 2,400)*Hz2rps; 

% Find the frequency response
H=freqresp(sys,wf);  % This ends up being a 3 dimensional array!
H=reshape(H,size(wf));	% this makes it the same size as wf

 
% Plot the amplitude as a function of frequency in Hz

loglog(wf*rps2Hz, Omega1^2*abs(H)); grid on;
xlabel('Frequency, omega (Hz)')
ylabel('Amplitude, X1(omega)/x_st ')
title( 'Optimized Amplitude for den Hartog oscillator' );
legend(['Omega =', num2str(Omega1*rps2Hz), '  omega 2 =', num2str(omega2*rps2Hz), ...
      '    mass ratio =', num2str(mu)]);
hold on;

axis([5 15  1   50]);

end;