disp(' ');
disp(' ss_ss_beam_frf.m  ver 1.3  December 10, 2013');
disp(' by Tom Irvine  Email: tom@vibrationdata.com ');
disp(' ');
disp(' Normal Modes & Frequency Response Function');
disp(' ');
disp(' where the force is uniform across the beam length. ');
disp(' ');
%
close all;
%
clear A;
clear root;
clear f;
clear H;
clear HA;
clear HM;
clear HM2;
clear beta;
clear fn;
clear YY;
clear beta;
clear xc;
clear transfer;
clear ptrans;
clear fn;
clear part;
clear emm;
%
clear HM_disp;
clear HM_bending_moment;
clear HM_shear;
%
fig_num=1;
%
disp(' Enter number of modes ');
n = input(' ');
%
iu=1;
%
[ccc,A,I,L]=geometry_entry(iu);
%
disp(' ');
%
[E,rho,mu]=materials(iu);
%
rho=rho*A;   % mass per unit length
%
inertia=I;
EI=E*inertia;
EI_term = sqrt(EI/rho);
%
mass=rho*L;
%
sq_mass=sqrt(mass);
%
beta=zeros(n,1);
omegan=zeros(n,1);
fn=zeros(n,1);
part=zeros(n,1);
%
for i=1:n
    beta(i)=i*pi/L;
    omegan(i)=(beta(i))^2*EI_term;
    fn(i)=omegan(i)/(2*pi);
end
%
for i=1:n
    part(i)=-sqrt(2*rho*L)*( cos(i*pi) - 1  )/(i*pi);
end
%
emm=part.^2;
emm=emm*386;
%
clear x;
clear y;
dx=L/200;
ax=0;
AC=sqrt(2)/sqrt(rho*L);
%
mode_num=1;
%
nu=201;
x=zeros(nu,1);
y=zeros(nu,1);
%
for i=1:nu
    x(i)=(i-1)*dx;
    y(i)=AC*sin(mode_num*pi*x(i)/L);
    ax=ax+(y(i))^2;
end
%
figure(fig_num);
fig_num=fig_num+1;
%
plot(x,y);
title('Mode 1');
xlabel('x (inch) ');
grid on;
%
disp(' ');
disp(' mode   fn(Hz)        PF     Effective Modal Mass (lbm)');
for i=1:n
    out1=sprintf('   %d %8.4g Hz\t %8.4g\t %8.4g ',i,fn(i),part(i),emm(i));
    disp(out1);
end
%
disp(' ');
disp(' Total Modal Mass ')
out1=sprintf('  %8.4g lbm ',sum(emm));
disp(out1);
%
[damp]=damping_entry(n);
%
nf=20000;
f(1)=1;
for k=2:nf
    f(k)=f(k-1)*2^(1/48);   % ignore warning
    if(f(k)>2000)
        break;
    end    
end
nf=k;
%
disp(' ');
clear x;
x=input(' Enter the response location (inch) ');
%
clear YY;
clear YYdd;
%
% i = mode_num
%
[YY,YYd,YYdd,YYddd]=ss_beam_modes(n,x,L,AC);
%
[HM_disp,HM_bending_moment,HM_shear_force]=...
                        beam_frf(n,nf,f,fn,damp,part,rho,EI,YY,YYdd,YYddd);
%
HM_bending_stress=abs(HM_bending_moment)*(ccc/I);
HM_average_shear_stress=abs(HM_shear_force)/A;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
f=fix_size(f);
HM_disp=fix_size(HM_disp);
HM_bending_moment=fix_size(HM_bending_moment);
HM_bending_stress=fix_size(HM_bending_stress);
HM_shear_force=fix_size(HM_shear_force);
HM_average_shear_stress=fix_size(HM_average_shear_stress);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(fig_num);
fig_num=fig_num+1;
plot(f,HM_disp);  
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
out1=sprintf('Displacement Frequency Response Function at x=%8.4g inch',x);
title(out1);
ylabel('[Displacement/(Force/Length)] [in/(lbf/in)] ');
xlabel('Frequency (Hz)');
grid on;
disp(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(fig_num);
fig_num=fig_num+1;
plot(f,HM_bending_moment);  
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
out1=sprintf('Bending Moment Frequency Response Function at x=%8.4g inch',x);
title(out1);
ylabel('[Bending Moment/(Force/Length)] [in-lbf/(lbf/in)] ');
xlabel('Frequency (Hz)');
grid on;
disp(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(fig_num);
fig_num=fig_num+1;
plot(f,HM_bending_stress);  
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
out1=sprintf('Bending Stress Frequency Response Function at x=%8.4g inch',x);
title(out1);
ylabel('[Stress/(Force/Length)] [psi/(lbf/in)] ');
xlabel('Frequency (Hz)');
grid on;
disp(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(fig_num);
fig_num=fig_num+1;
plot(f,HM_shear_force);  
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
out1=sprintf('Shear Force  Frequency Response Function at x=%8.4g inch',x);
title(out1);
ylabel('[Shear/(Force/Length)] [lbf/(lbf/in)] ');
xlabel('Frequency (Hz)');
grid on;
disp(' ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
out1=sprintf('\n Maximum values at %8.4g inch \n',x);
disp(out1);
%
xc=[f HM_disp];
disp_transfer=xc;
HM2_disp=HM_disp.*HM_disp;
clear C;
clear I;
[C,I]=max(xc);
xmax0=xc(I(2),2);
fmax0=xc(I(2),1);
%
ptrans=[f HM2_disp];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
xc=[f HM_bending_moment];
bending_moment_transfer=xc;        
clear C;
clear I;
[C,I]=max(xc);
xmax1=xc(I(2),2);
fmax1=xc(I(2),1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
xc=[f HM_bending_stress];
bending_stress_transfer=xc;        
clear C;
clear I;
[C,I]=max(xc);
xmax2=xc(I(2),2);
fmax2=xc(I(2),1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
xc=[f HM_shear_force];
shear_force_transfer=xc;        
clear C;
clear I;
[C,I]=max(xc);
xmax3=xc(I(2),2);
fmax3=xc(I(2),1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
xc=[f HM_average_shear_stress];
average_shear_transfer=xc;        
clear C;
clear I;
[C,I]=max(xc);
xmax4=xc(I(2),2);
fmax4=xc(I(2),1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
out0=sprintf('        Displacement:  max = %8.4g [in/(lbf/in)] at %8.4g Hz ',xmax0,fmax0);
out1=sprintf('      Bending Moment:  max = %8.4g [in-lbf/(lbf/in)] at %8.4g Hz',xmax1,fmax1);
out2=sprintf('      Bending Stress:  max = %8.4g [psi/(lbf/in)] at %8.4g Hz',xmax2,fmax2);
out3=sprintf('         Shear Force:  max = %8.4g [lbf/(lbf/in)] at %8.4g Hz',xmax3,fmax3);
out4=sprintf('Average Shear Stress:  max = %8.4g [psi/(lbf/in)] at %8.4g Hz',xmax4,fmax4);
%
disp(out0);
disp(out1);
disp(out2);
disp(out3);
disp(out4);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Output Arrays:');
disp(' ');
disp('         Displacement:  disp_transfer ');
disp('       Bending Moment:  bending_moment_transfer ');
disp('       Bending Stress:  bending_stress_transfer ');
disp('          Shear Force:  shear_force_transfer ');
disp(' Average Shear Stress:  average_shear_stress_transfer ');