%
disp(' ');
disp(' cant_beam_force_frf.m   ver 1.4   April 3, 2013 ');
disp(' ');
disp(' by Tom Irvine  Email: tomirvine@aol.com');
disp(' ');
disp(' Transverse vibration of a fixed-free beam, ');
disp(' subjected to force at its free end.     ');
disp(' Frequency response function output. ');
disp(' ');
%
close all;
%
clear beta;
clear bL;
clear omegan;
clear D;
clear H;
clear HM;
clear HVM;
clear HAM;
clear HL;
clear HLP;
clear HZM;
clear HZMP;
clear freq;
clear damp;
clear H_moment;
clear HM_moment;
clear HM_stress;
clear cc;
clear df;
%
tpi=2*pi;
%
nmodes=4;
%
[E,I,rho,c,L,A]=geometry_materials();
%
disp(' ');
cc=input(' Enter distance (in) from neutral axis to maximum stress locaton ');
%
disp(' ');
disp(' Enter the frequency step (Hz) ');
df=input(' ');
%
disp(' ');
disp(' Enter the maximum frequency (Hz) ');
fmax=input(' ');
%
x=L;
%
damp=input(' Enter the modal damping ');
%
disp(' ');
disp(' Enter the number of modes to include ');
nmodes=input(' ');
%
QM=1/(rho*A*L);
%
j=sqrt(-1);
%
beta(1)=1.87510;
beta(2)=4.69409;
beta(3)=5*pi/2;
beta(4)=7*pi/2;
%
for i=5:nmodes
    beta(i)=(2*i-1)*pi/2;
end
%
  beta=beta/L;
    bL=beta*L;
%
omegan=beta.*beta;
%
omegan=omegan*sqrt(E*I/(rho*A));
%
clear a;
for i=1:nmodes
    a=bL(i);
    D(i)=( cos(a)+cosh(a))/(sin(a)+sinh(a) );
end
%
disp(' ');
disp(' Natural Frequency (Hz) ');
for k=1:nmodes 
    out1=sprintf(' %8.4g ',omegan(k)/(2*pi));
    disp(out1);
end
%
%%%%%%%%%%%%%%%%%%%%%%%%
  Y = inline('(cosh(arg)-cos(arg))-Di*(sinh(arg)-sin(arg))','Di','arg')
DDY = inline('(cosh(arg)+cos(arg))-Di*(sinh(arg)+sin(arg))','Di','arg')
%%%%%%%%%%%%%%%%%%%%%%%%
%
clear freq;
clear H;
clear Hm;
clear Hp;
clear Hstrain;
clear Hstrain_mag;
clear HL;
%
disp(' ');
disp(' Calculating.... ');
%
for jk=1:100000
%
    f=jk*df;
%
    if(f>fmax)
        break;
    end
%
    freq(jk)=f;
    omega=tpi*f;
    H(jk)=0;
    H_moment(jk)=0;
%    
    for k=1:nmodes
%    
        omegad=omegan(k)*sqrt(1-damp^2);
        den=omegan(k)^2-omega^2 + j*2*damp*omega*omegan(k);
        arg1=bL(k);    
%
        x=L;
        arg2=beta(k)*x;
        Di=D(k);
        num=Y(Di,arg1)*Y(Di,arg2);
%
        H(jk)=H(jk)+ num/den;       
%
        x=0;
        arg2=beta(k)*x;
        num=Y(Di,arg1)*(beta(k))^2*DDY(Di,arg2);
        H_moment(jk)=H_moment(jk)+ num/den;   
%
    end
end
%
H=H*QM;
HM=abs(H);
%
H_moment=H_moment*QM*E*I;
HM_moment=abs(H_moment);
%
%%%%%%%%%%%%%%%
%
sz=size(freq);
if(sz(2)>sz(1))
    freq=freq';
end
%
sz=size(HM);
if(sz(2)>sz(1))
    HM=HM';
end
%
n=max(size(HM));
HVM=zeros(n,1);
HAM=zeros(n,1);
for i=1:n
    HVM(i)=tpi*freq(i)*HM(i);
    HAM(i)=tpi*freq(i)*HVM(i);
end
%
HAM=HAM/386;
%
sz=size(HM_moment);
if(sz(2)>sz(1))
    HM_moment=HM_moment';
end
%
HM_stress=HM_moment*cc/I;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(1);
plot(freq,HM);
xlabel(' Frequency (Hz) ');
ylabel(' Disp/Force (in/lbf) ');
title('Receptance Transfer Magnitude  x=L');
grid on;
%
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
% axis([ffmin,ffmax,ymin,ymax]);
%
figure(2);
plot(freq,HVM);
xlabel(' Frequency (Hz) ');
ylabel(' Vel/Force (in/sec/lbf) ');
title('Mobility Transfer Magnitude  x=L');
grid on;
%
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
%
figure(3);
plot(freq,HAM);
xlabel(' Frequency (Hz) ');
ylabel(' Accel/Force (G/lbf) ');
title('Accelerance Transfer Magnitude  x=L');
grid on;
%
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
%
figure(4);
plot(freq,HM_moment);
title(' Bending Moment Transfer Magnitude  x=0');
xlabel(' Frequency (Hz) ');
ylabel(' Moment/Force (in-lbf/lbf) ');
grid on;
%
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
% axis([ffmin,ffmax,ymin,ymax]);
%
figure(5);
plot(freq,HM_stress);
title(' Normal Stress Transfer Magnitude  x=0');
xlabel(' Frequency (Hz) ');
ylabel(' Stress/Force (psi/lbf) ');
grid on;
%
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
% axis([ffmin,ffmax,ymin,ymax]);
disp(' ');
disp(' Maximum values ');
disp(' ');
out1=sprintf(' displacement/force = %8.4g in/lbf      at x=L   ',max(HM));
out2=sprintf('     velocity/force = %8.4g in/sec/lbf  at x=L   ',max(HVM));
out3=sprintf(' acceleration/force = %8.4g G/lbf       at x=L   ',max(HAM));
%
out4=sprintf(' moment/force       = %8.4g in-lbf/lbf at x=0   ',max(HM_moment));
out5=sprintf(' stress/force       = %8.4g psi/lbf    at x=0   ',max(HM_stress));
%
disp(out1);
disp(out2);
disp(out3);
disp(' ');
disp(out4);
disp(out5);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear receptance;
clear mobility;
clear accelerance;
clear bending_moment;
clear normal_stress;
%
receptance=[freq HM];
mobility=[freq HVM];
accelerance=[freq HAM];
bending_moment=[freq HM_moment];
normal_stress=[freq HM_stress];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear receptance_power;
clear mobility_power;
clear accelerance_power;
clear bending_moment_power;
clear normal_stress_power;
%
receptance_power=[freq HM.*HM];
mobility_power=[freq HVM.*HVM];
accelerance_power=[freq HAM.*HAM];
bending_moment_power=[freq HM_moment.*HM_moment];
normal_stress_power=[freq HM_stress.*HM_stress];
%
while(1)
%
    disp(' ');
    disp(' Calculate response for a particular frequency case?');
    disp('  1=yes  2=no ');
    ic=input(' ');
%
    if(ic==1)
%
      disp(' ');
      freq=input(' Enter excitation frequency (Hz) ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
      clear omega;
      clear H;
      clear HVM;
      clear HAM;
      clear H_moment;
      clear HM_stress;
%
      omega=tpi*freq;
      H=0;
      H_moment=0;
%    
      for k=1:nmodes
%    
        omegad=omegan(k)*sqrt(1-damp^2);
        den=omegan(k)^2-omega^2 + j*2*damp*omega*omegan(k);
        arg1=bL(k);    
%
        x=L;
        arg2=beta(k)*x;
        Di=D(k);
        num=Y(Di,arg1)*Y(Di,arg2);
%
        H=H+ num/den;       
%
        x=0;
        arg2=beta(k)*x;
        num=Y(Di,arg1)*(beta(k))^2*DDY(Di,arg2);
        H_moment=H_moment+ num/den;   
%
      end
%
      H=H*QM;
      HM=abs(H);
%
      H_moment=H_moment*QM*E*I;
      HM_moment=abs(H_moment);
%
      HVM=tpi*freq*HM;
      HAM=tpi*freq*HVM;
%
      HAM=HAM/386;
%
      HM_stress=HM_moment*cc/I;   
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
     out1=sprintf(' displacement/force = %8.4g in/lbf      at x=L   ',max(HM));
     out2=sprintf('     velocity/force = %8.4g in/sec/lbf  at x=L   ',max(HVM));
     out3=sprintf(' acceleration/force = %8.4g G/lbf       at x=L   ',max(HAM));
%
     out4=sprintf(' moment/force       = %8.4g in-lbf/lbf at x=0   ',max(HM_moment));
     out5=sprintf(' stress/force       = %8.4g psi/lbf    at x=0   ',max(HM_stress));
%
     disp(' ');
     disp(out1);
     disp(out2);
     disp(out3);
     disp(' ');
     disp(out4);
     disp(out5);
%
    end
%
    if(ic==2)
        break;
    end
end