%
%  continuous_beam_base_frf.m  ver 1.5  January 29, 2013
%
function[fig_num,absolute_accel,relative_disp,Hmom]=...
   continuous_beam_base_frf(ModeShape,ModeShape_dd,sq_mass,beta,C,part,....
                                                fn,fig_num,x,damp,L,E,I,iu)
%
   clear f;
   nf=1000;
   n=length(fn);
%
   disp(' ');
   fmax=input(' Enter max freq (Hz) for frf plots:  ');
%
   f=zeros(nf,1);
   f(1)=1;
   for k=2:nf
        f(k)=f(k-1)*2^(1/48);
        if(f(k)>fmax)
            break;
        end    
   end
   nf=max(size(f));
   YY=zeros(n,1);
   for i=1:n
        arg=beta(i)*x;
        YY(i)=ModeShape(arg,C(i),sq_mass);
   end
%
   H=zeros(nf,1);
   H_relA=zeros(nf,1);
   for k=1:nf 
        H(k)=0;
        HrelA(k)=0;
        om=2*pi*f(k);
        for i=1:n
            pY=part(i)*YY(i);
            omn=2*pi*fn(i);
            num=-pY;
            den=(omn^2-om^2)+(1i)*2*damp(i)*om*omn;
            H(k)=H(k)+num/den;
%
        end
        H_relA(k)=-om^2*H(k);
   end
%
    H=abs(H);
    HA=abs(H_relA+1);   
%
    if(iu==1)
        H=H*386;
    else
        H=H*9.81;       
    end
%
    maxH=0.;
    maxHA=0.;
    maxF=0.;
    maxFA=0.;
%
    out1=sprintf('\n nf=%d \n',nf);
    disp(out1);
%
    for k=1:nf
        if(H(k)>maxH)
            maxH=H(k);
            maxF=f(k);
        end
        if(HA(k)>maxHA)
            maxHA=HA(k);
            maxFA=f(k);
        end
    end       
%
    sz=size(f);
    if(sz(2)>sz(1))
        f=f';
    end
    sz=size(H);
    if(sz(2)>sz(1))
        H=H';
    end 
    sz=size(HA);
    if(sz(2)>sz(1))
        HA=HA';
    end 
%
    n=length(f);
    for i=n:-1:1
%
       if(f(i)==0)
            f(i)=[];
            H(i)=[];
           HA(i)=[];
       H_relA(i)=[];  
       end
%
    end
%
    figure(fig_num);
    fig_num=fig_num+1;
    plot(f,H);
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
%
    if(iu==1)
        out1=sprintf('Frequency Response Function  x=%7.4g inch',x);
        ylabel('Rel Disp (in)/ Base Accel (G)');
    else
        out1=sprintf('Frequency Response Function  x=%7.4g m',x);
        ylabel('Rel Disp (m)/ Base Accel (G)');        
    end
%
%
    title(out1);
    xlabel('Frequency (Hz)');
    set(gca,'FontSize',12);
    h = get(gca, 'xlabel');
    set(h, 'FontName', 'Arial','FontSize',12);
    h = get(gca, 'ylabel');
    set(h, 'FontName', 'Arial','FontSize',12); 
    h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',12);
    grid on;
    disp(' ');
    if(iu==1)
        out1=sprintf('  max Rel Disp FRF = %9.5f (in/G) at %8.4g Hz ',maxH,maxF);
    else
        out1=sprintf('  max Rel Disp FRF = %9.5f (m/G) at %8.4g Hz ',maxH,maxF);       
    end
    disp(out1);
%
    clear relative_disp;
    relative_disp=[f H];
%
    figure(fig_num);
    fig_num=fig_num+1;
    plot(f,abs(H_relA));
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
    if(iu==1)
        out1=sprintf('Frequency Response Function  x=%7.4g inch',x);
    else
        out1=sprintf('Frequency Response Function  x=%7.4g m',x);        
    end
    title(out1);
    ylabel('Relative Accel (G)/ Base Accel (G)');
    xlabel('Frequency (Hz)');
    set(gca,'FontSize',12);
    h = get(gca, 'xlabel');
    set(h, 'FontName', 'Arial','FontSize',12);
    h = get(gca, 'ylabel');
    set(h, 'FontName', 'Arial','FontSize',12); 
    h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',12);    
    grid on;
%
    figure(fig_num);
    fig_num=fig_num+1;
    plot(f,HA);
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
%    
    if(iu==1)
        out1=sprintf('Frequency Response Function  x=%7.4g inch',x);
    else
        out1=sprintf('Frequency Response Function  x=%7.4g m',x);        
    end
%    
    title(out1);
    ymax= 10^(ceil(log10(max(abs(HA)))+0.1));
    ymin= 10^(floor(log10(min(abs(HA)))-0.1));
%
    LL=ymax/10000;
%
    if(ymin<LL)
        ymin=LL;
    end
%
    axis([min(xlim),max(xlim),ymin,ymax]);
    ylabel('Absolute Accel / Base Accel ');
    xlabel('Frequency (Hz)');
    set(gca,'FontSize',12);
    h = get(gca, 'xlabel');
    set(h, 'FontName', 'Arial','FontSize',12);
    h = get(gca, 'ylabel');
    set(h, 'FontName', 'Arial','FontSize',12); 
    h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',12);
    grid on;
    disp(' ');
    out1=sprintf('  max Absolute Accel FRF = %8.4g (G/G) at %8.4g Hz ',maxHA,maxFA);
    disp(out1);
%
    clear absolute_accel;
    absolute_accel=[f HA];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear x;
x=[0 0.5*L L];
%
nx=length(x);
%
n=length(fn);
%
nf=length(f);
%
conv=E*I;
if(iu==1)
    conv=conv*386;
else
    conv=conv*9.81;
end    
%
HMM=zeros(nf,nx);
%
for j=1:nx
%
    for i=1:n 
%
        arg=beta(i)*x(j);
        num=-part(i)*ModeShape_dd(arg,C(i),beta(i),sq_mass);
%
        omn=2*pi*fn(i);
%
        for k=1:nf
            om=2*pi*f(k);
%
            den=(omn^2-om^2)+(1i)*2*damp(i)*om*omn;
            HMM(k,j)=HMM(k,j)+(conv*num/den);
%
        end      
%
    end
%
end
%
HMM=abs(HMM);
%
Hmom=[f HMM];
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
for j=1:3
%
   if(max(HMM(:,j))>1.0e-10)
        figure(fig_num);
        fig_num=fig_num+1;
        plot(f,HMM(:,j));
        set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
        if(iu==1)
            out1=sprintf('Frequency Response Function  x=%7.4g inch',x(j));
        else
            out1=sprintf('Frequency Response Function  x=%7.4g m',x(j));            
        end
        title(out1);
        ymax= 10^(ceil(log10(max(abs(HMM(:,j))))+0.1));
        ymin= 10^(floor(log10(min(abs(HMM(:,j))))-0.1));
%
        LL=ymax/10000;
%
        if(ymin<LL)
            ymin=LL;
        end
%
        axis([min(xlim),max(xlim),ymin,ymax]);
        if(iu==1)
            ylabel('Bending Moment (in lbf)/ Base Accel(G) ');
        else
            ylabel('Bending Moment (N m/ Base Accel(G) ');            
        end
        xlabel('Frequency (Hz)');
        set(gca,'FontSize',12);
        h = get(gca, 'xlabel');
        set(h, 'FontName', 'Arial','FontSize',12);
        h = get(gca, 'ylabel');
        set(h, 'FontName', 'Arial','FontSize',12); 
        h = get(gca, 'title');
        set(h, 'FontName', 'Arial','FontSize',12);
        grid on;
        disp(' ');     
   end
%
end