%
disp(' ');
disp(' three_parameter_isolation.m  ver 1.0  November 29, 2012 ');
disp(' ');
disp(' by Tom Irvine ');
disp(' ');
%
close all;
clear f;
clear A;
clear A_mag;
clear A_phase;
%
tpi=2*pi;
%
 m=input(' Enter  m: ');
 c=input(' Enter  c: ');
k1=input(' Enter k1: ');
k2=input(' Enter k2: ');
%
omega_1=sqrt(k1/m);
omega_2=sqrt(k2/m);
d=c/m;
%
minf=(0.2)*min([omega_1 omega_2])/tpi;
minf=10^(floor(log10(minf)));
%
maxf=5*max([omega_1 omega_2])/tpi;
maxf=10^(ceil(log10(maxf)));
%
freq=minf;
for i=1:1000
    omega=tpi*freq;
    s=(1i)*omega;
%
    term= (  omega_2^2 + omega_1^2 )*d*s + omega_1^2*omega_2^2;
    num=term;
    den= d*s^3 + omega_2^2*s^2 + term;
    A(i)=num/den;
    f(i)=freq;
%    
    if(freq==maxf)
        break;
    end
    freq=freq*2^(1/48);
    if(freq>maxf)
        freq=maxf;
    end
end
%
A_mag=abs(A);
A_phase=atan2(imag(A),real(A))*180/pi;
%
C=[f' A_mag'];  
%
[zmax,fz_max]=find_max(C);
%
out5 = sprintf('\n Peak is %10.5g (G/G) at %10.5g Hz ',zmax,fz_max);
disp(out5);
%
figure(1)
%
subplot(3,1,1);
plot(f,A_phase);
title('Transmissibility');
xlabel(' Frequency (Hz)');
ylabel(' Phase (deg)'); 
grid on; 
axis([f(1),maxf,-180,0]); 
set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale',...
                         'log','YScale','lin','ytick',[-180,-90,0]);
%
subplot(3,1,[2 3]);
plot(f,A_mag);
ylabel('Trans(G/G)');
xlabel('Frequency(Hz)');
grid on;
set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log'); 