% resonator frf computations  
% by k c park  10 april 2005
% compute x_2/x_1st;  x1st = F1/K1;
clear all;

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

% Model parameters
K1 = 2720;
K12 = 75.868;
M1 = 5.73266e-13;

omega_b =  Hz2rps*(10.475e+6);
Mb = 5.72e-13;



zeta =[0  0.001  0.005  0.01  0.05];

for i =1:length(zeta),  %zeta is the damping ratio
%
% Define a system in transfer function form
%
C = 2*zeta(i)*Mb*sqrt(K1/Mb);
%
a0=0; a1 =0;  a2 = 0; a3=0; a4 = K1*K12;
%
b0 = M1^2; 
b1 = 2*M1*C; 
b2 = C^2 + 2*(K1+K12)*M1;
b3 = 2*(K1+K12)*C; 
b4 = (K1+K12)^2 - K12^2;

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

% Pick a range of frequencies to look at from 0 Hz to 100 Hz
wf= logspace( 6.95, 7.15, 10^4)*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

semilogy(wf*rps2Hz, abs(H)); grid on;

xlabel(' Driving Frequency  (Hz)')
ylabel('Amplitude, X2(omega)/x1_st ')
title( 'Double Beam Resonator' );
%legend(['Omega =', num2str(Omega1*rps2Hz), '  omega 2 =', num2str(omega2*rps2Hz), ...
%     '    mass ratio =', num2str(mu)]);
hold on;

% this inverval draws frequency from .8*10^7  to 1.5*10^7 range
axis([0.9*10^7   1.3*10^7  1.e-1   1.e+6]);
%
%   For 0.9*10^7 to 10^7 interval use this:
%  axis([1.07*10^7  1.14*10^7   1.e-1   1.e+6]);
%
end;