disp(' ');
disp(' fixed_fixed_beam_spatial_frf.m  ver 1.1  July 2, 2012');
disp(' by Tom Irvine  Email: tomirvine@aol.com ');
disp(' ');
disp(' Normal Modes & Frequency Response Function');
disp(' ');
disp(' where the force distribution is equal to the fundamental mode. ');
disp(' ');
%
close all;
%
clear root;
clear f;
clear ccc;
clear H;
clear HA;
clear HM;
clear HM2;
clear beta;
clear YY;
clear beta;
clear xc;
clear transfer;
clear ptrans;
clear C;
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
%
clear EI;
clear EI_term;
inertia=I;
EI=E*inertia;
EI_term = sqrt(EI/rho);
%
mass=rho*L;
%
sq_mass=sqrt(mass);
%
root(1)=4.73004;
root(2)=7.8532;
root(3)=10.9956;
root(4)=14.13717;
root(5)=17.27876;
%
for(i=1:n)
    if(i>=6)
        root(i)=(0.5 + i)*pi;
    end  
    beta(i)=root(i)/L;
    omegan(i)=(beta(i))^2*EI_term;
    fn(i)=omegan(i)/(2*pi);
end
%
for(i=1:n)
    bL=root(i);
    C(i)=(sinh(bL)+sin(bL))/(cosh(bL)-cos(bL));
end
%
for(i=1:n)
    arg=root(i);
    p2=(sinh(arg)-sin(arg))-C(i)*(cosh(arg)+cos(arg));
    arg=0;
    p1=(sinh(arg)-sin(arg))-C(i)*(cosh(arg)+cos(arg));
    part(i)=(p2-p1)/beta(i);
end
part=part/L;
emm=part.^2;
emm=emm*mass*386;
part=part*sq_mass;
%
clear x;
clear y;
clear yydd;
dx=L/200;
ax=0;
for(i=1:201)
    x(i)=(i-1)*dx;
    arg=beta(1)*x(i);
    y(i)=                (cosh(arg)-cos(arg))-C(1)*(sinh(arg)-sin(arg));
    yydd(i)=(beta(1))^2*((cosh(arg)+cos(arg))-C(1)*(sinh(arg)+sin(arg))); 
end
y=      y/sqrt(rho*L);
yydd=yydd/sqrt(rho*L);
%
figure(fig_num);
fig_num=fig_num+1;
plot(x,y);
title('Mode 1');
xlabel('x (inch) ');
grid on;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(fig_num);
fig_num=fig_num+1;
plot(x,yydd);
title('Second derivative of 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(' ');
   x=input(' Enter the response location (inch) ');
%
   clear YY;
   clear YYdd;
   for(i=1:n)
        arg=beta(i)*x;
        YY(i)  =             (cosh(arg)-cos(arg))-C(i)*(sinh(arg)-sin(arg));
        YYdd(i)=(beta(i))^2*((cosh(arg)+cos(arg))-C(i)*(sinh(arg)+sin(arg)));   
   end
%
     YY=  YY/sq_mass;
   YYdd=YYdd/sq_mass;
%
    clear HM;
    clear HM_stress;
%
    n=1;
%
    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; 
            arg=beta(i)*x;
%            
            num=(cosh(arg)-cos(arg))-C(i)*(sinh(arg)-sin(arg));
            HM(k)=HM(k)+num/den;
%
            num=(beta(i))^2*((cosh(arg)+cos(arg))-C(i)*(sinh(arg)+sin(arg)));      
            HM_stress(k)=HM_stress(k)+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 ');