disp('    sdof_dof_force_frf.m   ver 1.5  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 an applied force. ');
%
close all hidden;
%
clear k;
clear m;
clear Eigenvalues;
clear ModeShapes;
clear omega;
clear omegan;
clear fn;
clear MS;
clear MST;
clear ModalMass;
clear part;
clear QTMQ;
clear r;
clear H;
clear H1;
clear H2;
%
clear acc_m;
clear mob_m;
clear recp_m;
%
clear accelerance;
clear mobility;
clear receptance;
%
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';
    receptance_unit='in/lbf';
    mobility_unit='in/sec/lbf';
    accelerance_unit='G/lbf';
else
    mass_unit='kg';
    stiffness_unit='N/m';
    force_unit='N';
    receptance_unit='m/N';
    mobility_unit='m/sec/N';
    accelerance_unit='G/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(' ');
disp('  Natural Frequency = ');
out1=sprintf('\n  %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(' ');
%
%  Find Receptance  H = X/F
%
clear f;
clear H;
clear Hm;
clear Ha;
%
num=round(fmax/df);
%
f=linspace(0,(num-1)*df,num);
%
f(1)=[];
num=num-1;
%
omega=tpi*f;
%
H=zeros(num,1);
Hm=zeros(num,1);
Ha=zeros(num,1);
%
for i=1:num
     H(i)=(omegan^2/k)/((omegan^2-omega(i)^2)+(1i)*(2*damp*omega(i)*omegan));
    Hm(i)=(1i)*omega(i)*H(i);
    Ha(i)=(1i)*omega(i)*Hm(i);    
%
end
%
%
if(iu==1)
    Ha=Ha/386;    
end
%
omegad=omegan*sqrt(1-damp^2);
pole1=-damp*omegan+(1i)*omegad;
pole2=-damp*omegan-(1i)*omegad;
disp(' ');
out1=sprintf('pole1:  %8.4g + j%8.4g  rad/sec',real(pole1),imag(pole1));
out2=sprintf('pole2:  %8.4g - j%8.4g  rad/sec',real(pole2),abs(imag(pole2)));
disp(out1);
disp(out2);
disp(' ');
%
omega=omegan;
clear Hfn;
Hfn=(omegan^2/k)/((omegan^2-omega^2)+(1i)*(2*damp*omega*omegan));
%
figure(1);
%
plot(f,abs(H));
grid;
set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log');
%
title(' Receptance FRF (disp/force)');
xlabel('Frequency (Hz)');
out1=sprintf(' Receptance (%s)',receptance_unit);
ylabel(out1);
    h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',11)
        h = get(gca, 'xlabel');
    set(h, 'FontName', 'Arial','FontSize',11)
         h = get(gca, 'ylabel');
    set(h, 'FontName', 'Arial','FontSize',11) 
%
%
figure(2);
%
plot(f,abs(Hm));
grid;
set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log');
%
title(' Mobility FRF (velocity/force)');
xlabel('Frequency (Hz)');
out1=sprintf(' Mobility (%s)',mobility_unit);
ylabel(out1);
    h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',11)
        h = get(gca, 'xlabel');
    set(h, 'FontName', 'Arial','FontSize',11)
         h = get(gca, 'ylabel');
    set(h, 'FontName', 'Arial','FontSize',11) 
%
figure(3);
%
plot(f,abs(Ha));
grid;
set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log');
%
title(' Accelerance FRF (accel/force)');
xlabel('Frequency (Hz)');
%
out1=sprintf(' Accelerance (%s)',accelerance_unit);
ylabel(out1);
%
    h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',11)
        h = get(gca, 'xlabel');
    set(h, 'FontName', 'Arial','FontSize',11)
         h = get(gca, 'ylabel');
    set(h, 'FontName', 'Arial','FontSize',11) 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(4);
%
plot(f,abs(H));
grid;
set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log');
%
a=1;
b=2*omegan^2*(-1+2*damp^2);
c=omegan^4*(1-8*damp^2);
%
om2=(-b+sqrt(b^2-4*a*c))/(2*a);   %  Half-power bandwidth frequencies
om1=(-b-sqrt(b^2-4*a*c))/(2*a);
%
om2=sqrt(om2);
om1=sqrt(om1);
%
clear Hfn1;
Hfn1=(omegan^2/k)/((omegan^2-om1^2)+(1i)*(2*damp*om1*omegan));
clear Hfn2;
Hfn2=(omegan^2/k)/((omegan^2-om2^2)+(1i)*(2*damp*om2*omegan));
%
plot(real(H),imag(H));
hold on
plot(real(Hfn),imag(Hfn),'x','MarkerSize',12);
% plot(real(Hfn1),imag(Hfn1),'d','MarkerSize',12);
% plot(real(Hfn2),imag(Hfn2),'d','MarkerSize',12);
%
out1=sprintf('Nyquist FRF of Receptance  (%s)',receptance_unit);
title(out1);
%
ylabel('Imag');
xlabel('Real');
delta=0.2*max(real(H));
xmin=1*min(real(H));
xmax=1*max(real(H));
ymin=-(xmax-xmin);
ymax=0;
%
xmin=xmin-delta;
xmax=xmax+delta;
ymin=ymin-delta;
ymax=ymax+delta;
%
axis([xmin,xmax,ymin,ymax]);
axis square
grid;
hold off
    h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',11)
        h = get(gca, 'xlabel');
    set(h, 'FontName', 'Arial','FontSize',11)
         h = get(gca, 'ylabel');
    set(h, 'FontName', 'Arial','FontSize',11) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(5)
plot(real(Hm),imag(Hm));
hold on
%%% plot(real(Hfn),imag(Hfn),'x','MarkerSize',12);
% plot(real(Hfn1),imag(Hfn1),'d','MarkerSize',12);
% plot(real(Hfn2),imag(Hfn2),'d','MarkerSize',12);
%
%
out1=sprintf('Nyquist FRF of Mobility  (%s)',mobility_unit);
title(out1);
ylabel('Imag');
xlabel('Real');
delta=0.2*max(imag(Hm));
ymin=min(imag(Hm));
ymax=max(imag(Hm));
xmin=0.;
xmax=ymax-ymin;
%
xmin=xmin-delta;
xmax=xmax+delta;
ymin=ymin-delta;
ymax=ymax+delta;
%
axis([xmin,xmax,ymin,ymax]);
axis square
grid;
hold off
    h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',11)
        h = get(gca, 'xlabel');
    set(h, 'FontName', 'Arial','FontSize',11)
         h = get(gca, 'ylabel');
    set(h, 'FontName', 'Arial','FontSize',11) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(6)
receptance(:,1)=f;
receptance(:,2)=H;
mobility(:,1)=f;
mobility(:,2)=Hm;
accelerance(:,1)=f;
accelerance(:,2)=Ha;
%
%
plot(real(Ha),imag(Ha));
hold on
%%% plot(real(Hfn),imag(Hfn),'x','MarkerSize',12);
% plot(real(Hfn1),imag(Hfn1),'d','MarkerSize',12);
% plot(real(Hfn2),imag(Hfn2),'d','MarkerSize',12);
%
out1=sprintf('Magnitude Accelerance (%s)',accelerance_unit);
title(out1);
ylabel('Imag');
xlabel('Real');
%
delta=0.2*max(real(Ha));
xmin=1*min(real(Ha));
xmax=1*max(real(Ha));
ymin=0;
ymax=(xmax-xmin);
%
xmin=xmin-delta;
xmax=xmax+delta;
ymin=ymin-delta;
ymax=ymax+delta;
%
axis([xmin,xmax,ymin,ymax]);
axis square
grid;
hold off
    h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',11)
        h = get(gca, 'xlabel');
    set(h, 'FontName', 'Arial','FontSize',11)
         h = get(gca, 'ylabel');
    set(h, 'FontName', 'Arial','FontSize',11) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(7);
ph=unwrap(angle(H));
ph=ph*180/pi;
plot(ph,abs(H));
title(' Nichols Plot ');
xlabel('Phase [Degrees]');
out1=sprintf('Magnitude Receptance (%s)',receptance_unit);
ylabel(out1);
grid on;
%
    h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',11)
        h = get(gca, 'xlabel');
    set(h, 'FontName', 'Arial','FontSize',11)
         h = get(gca, 'ylabel');
    set(h, 'FontName', 'Arial','FontSize',11) 
%
%
sz=size(f);
if(sz(2)>sz(1))
    f=f';
end
%
acc_m =[f abs(accelerance(:,2))];
mob_m =[f abs(mobility(:,2))];
recp_m=[f abs(receptance(:,2))];
%
disp(' ');
disp('  Accelerance Magnitude ');
[ma,ff]=find_max(acc_m);
out1=sprintf('\n max=%8.4g %s at %8.4g Hz \n',ma,accelerance_unit,ff);
disp(out1);
%
disp('  Mobility Magnitude ');
[ma,ff]=find_max(mob_m);
out1=sprintf('\n max=%8.4g %s at %8.4g Hz \n',ma,mobility_unit,ff);
disp(out1);
%
disp('  Receptance Magnitude ');
[ma,ff]=find_max(recp_m);
out1=sprintf('\n max=%8.4g %s at %8.4g Hz \n',ma,receptance_unit,ff);
disp(out1);
%
%
disp(' ');
disp(' Output files ');
disp('   receptance ');
disp('   mobility ');
disp('   accelerance ');