disp('    sdof_dof_base_frf.m   ver 1.0  September 24, 2012 ');
disp(' ');
disp('    by Tom Irvine   Email: tomirvine@aol.com ');
disp(' ');
disp(' This program finds the natural frequency for a ');
disp(' single-degree-of-freedom system.');
disp(' The equation of motion is:   M (d^2x/dt^2) + K x = 0 ');
disp(' ');
disp(' It also finds the frequency response function for base excitation. ');
%
clear k;
clear m;
clear omega;
clear omegan;
clear fn;
clear MS;
clear MST;
clear ModalMass;
clear part;
clear H;
%
clear acc_m;
clear rd_m
%
clear acceleration_trans;
clear relative_disp_trans;
%
tpi=2.*pi;
%
disp(' Select units: 1=English  2=metric ');
iunit=input(' ');
iu=iunit;
%
if(iunit==1)
    mass_unit='lbm';
    stiffness_unit='lbf/in';
    force_unit='lbf';
else
    mass_unit='kg';
    stiffness_unit='N/m';
    force_unit='N';
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    
out1=sprintf(' Enter mass (%s)',mass_unit);
disp(out1);
m=input(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if(iunit==1)
    m=m/386;
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
out1=sprintf(' Enter stiffness (%s)',stiffness_unit);
disp(out1);
k=input(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Enter modal damping ratio ');
damp=input(' ');
%
omegan = sqrt(k/m);
fn=omegan/tpi;
%
disp(' ');
out1=sprintf('\n Natural Frequency =  %8.4g Hz',fn);
disp(out1);
%
fmax=20*fn;
df=fn/200;
dff=df;
%
disp(' ')
out1=sprintf(' Enter frequency step Hz (suggest < %6.3g )  ',df);
disp(out1);
df=input(' ');
%
clear f;
clear H;
clear H_acc;
clear H_rd;
clear H_tf;
%
num=round(fmax/df);
%
f=linspace(0,(num-1)*df,num);
%
f(1)=[];
num=num-1;
%
omega=tpi*f;
%
H_acc=zeros(num,1);
H_rd=zeros(num,1);
%
for i=1:num
%
    jterm=(1i)*2*damp*omega(i)*omegan;
%
    den=((omegan^2-omega(i)^2)+jterm);
%
    H_acc(i)=(omegan^2 +jterm )/den; 
    H_rd(i)= -1/den; 
%
end
%
sz=size(f);
if(sz(2)>sz(1))
    f=f';
end
%
if(iu==1)
    H_rd=H_rd*386;
end
%
figure(1);
plot(f, abs(H_acc));
title(' Acceleration Transfer Magnitude ');
if(iu==1)
    ylabel('Trans(G/G)');
else
    ylabel('Trans((m/sec^2)/(m/sec^2))');
end
grid on;
set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log');
%
figure(2);
plot(f, abs(H_rd));
title(' Relative Displacement Transfer Magnitude ');
if(iu==1)
    ylabel('Trans(in/G)');
else
    ylabel('Trans(m/(m/sec^2))');
end
xlabel('Frequency(Hz)');
grid on;
set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log');
%
acceleration_trans=[f  H_acc];
relative_disp_trans=[f H_rd];
%
acc_m=[f  abs(H_acc)];
rd_m=[f abs(H_rd)];
%
disp('  Acceleration Transfer Magnitude ');
[ma,ff]=find_max(acc_m);
out1=sprintf('\n max=%8.4g at %8.4g Hz \n',ma,ff);
disp(out1);
%
disp('  Relative Displacement Transfer Magnitude ');
[mrd,ff]=find_max(rd_m);
out1=sprintf('\n max=%8.4g at %8.4g Hz \n',mrd,ff);
disp(out1);
%
disp(' ');
disp(' Output Files: ');
if(iu==1)
    disp('     acceleration_trans (Hz & G/G)');    
    disp('  relative displacement (Hz & in/G)');
else
    disp('     acceleration_trans (Hz & (m/sec^2)/(m/sec^2)');    
    disp('  relative displacement (Hz & m/(m/sec^2)');    
end
%