disp(' ');
disp(' ss_ss_beam_spatial_frf.m  ver 1.2  July 2, 2012');
disp(' by Tom Irvine  Email: tom@vibrationdata.com ');
disp(' ');
disp(' Normal Modes & Frequency Response Function');
disp(' ');
disp(' where the force shape is the same as the normalize first mode ');
disp(' across the beam length.  i.e. half-sine shape');
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;
%
fig_num=1;
%
disp(' Enter number of modes ');
n = input(' ');
%
iu=1;
%
[ccc,A,I,L]=geometry_entry(iu);
%
[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);
%
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;
%
for(i=1:201)
    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);
%
    x=0;
    j=sqrt(-1);   
%
    nf=20000;
    f(1)=1;
    for(k=2:nf)
        f(k)=f(k-1)*2^(1/48);
        if(f(k)>2000)
            break;
        end    
    end
   nf=k;
%
   disp(' ');
   clear x;
   x=input(' Enter the response location (inch) ');
%
   clear YY;
%
% i = mode_num
%
    n=1;  % only fundamental mode active
%
    clear HM;
    clear HM_stress;
%
    for(k=1:nf)
        HM(k)=0;
        HM_stress(k)=0;
        om=2*pi*f(k);
        for(i=1:n)
            omn=2*pi*fn(i);
            den=(omn^2-om^2)+j*2*damp(i)*om*omn;
%            
            num=sin(pi*x/L);
            HM(k)=HM(k)+num/den;
 %
            HM_stress(k)=HM_stress(k)-(pi/L)^2*num/den;
 %
        end
    end
   %
%%
    HM=abs(HM)/rho;   %% divide by mass/length
    HM_stress=abs(HM_stress)*(EI/rho)*(ccc/I);
%
    figure(fig_num);
    fig_num=fig_num+1;
    plot(f,HM);  
    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_stress);  
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
    out1=sprintf('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(' ');
%
    f=f';
    HM=HM';
    xc=[f HM];
    disp_transfer=xc;
    HM2=HM.*HM;
    clear C;
    clear I;
    [C,I]=max(xc);
    xmax=xc(I(2),2);
    fmax=xc(I(2),1);
%
    out1=sprintf('Displacement:  max = %8.4g [in/(lbf/in)] at %8.4g Hz  ',xmax,fmax);
    disp(out1);
%
    ptrans=[f HM2];
%
    HM_stress=HM_stress';
    xc=[f HM_stress];
    stress_transfer=xc;        
    clear C;
    clear I;
    [C,I]=max(xc);
    xmax=xc(I(2),2);
    fmax=xc(I(2),1);
%
    disp(' ');
    out1=sprintf('Stress:  max = %8.4g [psi/(lbf/in)] at %8.4g Hz  ',xmax,fmax);
    disp(out1);
%
    disp(' ');
    disp(' The displacement output array name is:  disp_transfer ');
%
    disp(' ');
    disp(' The stress output array name is:  stress_transfer ');