%
disp(' ');
disp(' long_force_frf.m   ver 1.0   October 26, 2010 ');
disp(' ');
disp(' Longitudinal vibration of a fixed-free rod, ');
disp(' subjected to force.  Frequency response function output. ');
disp(' ');
%
[E,I,rho,c,L,A]=geometry_materials();
%
disp(' ');
disp(' Enter the maximum frequency (Hz) ');
fmax=input(' ');
%
x=L;
%
damp=input(' Enter the modal damping ');
%
QM=2/(rho*A*L);
%
j=sqrt(-1);
%
disp(' ');
disp(' Natural Frequency (Hz) ');
for(k=1:3)
%
    i=2*k-1;
%    
    omegan(k)=i*pi*c/(2*L);
    out1=sprintf(' %8.4g ',omegan(k)/(2*pi));
    disp(out1);
%
end
%
clear freq;
clear H;
clear Hm;
clear Hp;
clear Hstrain;
clear Hstrain_mag;
clear HL;
%
for(jk=1:10000)
    f=jk;
    if(f>fmax)
        break;
    end
    freq(jk)=f;
    omega=2*pi*f;
    H(jk)=0;
    Hstrain(jk)=0;
%    
    for(k=1:3)
%    
        i=2*k-1;
%    
        omegad=omegan(k)*sqrt(1-damp^2);
        den=omegan(k)^2-omega^2 + j*2*damp*omega*omegan(k);
        arg1=omegan(k)*L/c;    
%
        x=L;
        arg2=omegan(k)*x/c;
        num=sin(arg1)*sin(arg2);
%
        H(jk)=H(jk)+ num/den;
%
        x=0;
        arg2=omegan(k)*x/c;
        num=sin(arg1)*cos(arg2);
%
        Hstrain(jk)=Hstrain(jk)+ (omegan(k)/c)*(num/den);        
%
    end
end
%
H=H*QM;
HM=abs(H);
%
Hstrain=Hstrain*QM;
Hstrain_mag=abs(Hstrain);
%
figure(1);
plot(freq,HM);
xlabel(' Frequency (Hz) ');
ylabel(' H (in/lbf) ');
title('x=L');
grid on;
%
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
% axis([ffmin,ffmax,ymin,ymax]);
%%%%%%%%%%%%%%%
%
sz=size(freq);
if(sz(2)>sz(1))
    freq=freq';
end
%
sz=size(HM);
if(sz(2)>sz(1))
    HM=HM';
end
%
HL=[freq HM];
HLP=[freq HM.*HM];
%
%%%%%%%%%%%%%%%
%
figure(2);
plot(freq,Hstrain_mag);
xlabel(' Frequency (Hz) ');
ylabel(' H (strain/force) ');
title('x=0');
grid on;
%
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
% axis([ffmin,ffmax,ymin,ymax]);
%%%%%%%%%%%%%%%
%
sz=size(freq);
if(sz(2)>sz(1))
    freq=freq';
end
%
sz=size(Hstrain_mag);
if(sz(2)>sz(1))
    Hstrain_mag=Hstrain_mag';
end
%
HS=[freq Hstrain_mag];
HSP=[freq Hstrain_mag.*Hstrain_mag];
%