disp(' ')
disp(' circular_SS.m   ver 2.3  March 1, 2012 ')
disp(' ')
disp(' by Tom Irvine   Email: tomirvine@aol.com ')
disp(' ');
disp(' This program calculates the natural frequencies of the ');
disp(' bending modes of a simply-supported circular plate ');
disp(' The solution is derived from Bessel functions. ');
disp(' It also calculates the base excitation transfer function. ');
disp(' ');
%
close all hidden;
%
clear freq;
clear damp;
clear Hra;
clear Haa;
clear Hrd;
clear Haa_phase;
clear Hrd_phase;
clear Haa_mag;
clear Hrd_mag;
%
clear n;
clear fn;
clear a;
clear b;
clear q;
clear xo;
%
clear x;
clear y;
clear root;
clear rootn;
clear rootk;
%
clear PF;
%
clear Dc;
clear Cc;
%
tpi=2*pi;
%
[iu,D,rho,mass,radius,h,mu]=circular_plate_materials();
%
DD=sqrt(D/(rho*h));
%
delta=0.1;
%
num=6;
num2=num^2;
nm1=num-1;
%
for(ijk=1:num2)
    root(ijk)=0;
    rootn(ijk)=0;
    rootk(ijk)=0;
    Cc(ijk)=0;
    Dc(ijk)=0;
    PF(ijk)=0;
end
%
ijk=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
num_roots=(nm1+1)^2;
%
progressbar % Create figure and set starting time
%
for n=0:nm1
%
    clear a;
    clear b;
    clear x;
    clear y;
%
    j=1;
    delta=0.01;
    for(i=1:12000)
        xx=i*delta;
        x(i)=xx;
%        
        ya=Bessel_SS(n,xx,mu);   
%
        if(i>2)
            if(ya*yb<=0)
                b(j)=xx;
                j=j+1;
            end
        end
        yb=ya;    
        if(j>4)
            break;
        end
%        
    end
    a=b-delta;
%
%   secant
%
    k=0;
    for i=1:max(size(a))
%
        for j=1:150
%            
            xx=a(i);    
            ya=Bessel_SS(n,xx,mu);       
%
            xx=b(i);      
            yb=Bessel_SS(n,xx,mu);   
%
            if(abs(b(i)-a(i))<0.00001)
                break;
            end
%
            m=(yb-ya)/(b(i)-a(i));
            bb=ya-m*a(i);
            ab=-bb/m;
%
            xx=ab;           
            yab=Bessel_SS(n,xx,mu);   
%
            if(ya*yab<=0)
                b(i)=ab;    
            else
                a(i)=ab;
            end
%
        end
        root(ijk)=0.25*(a(i)+b(i))^2;
        rootn(ijk)=n;
        rootk(ijk)=k;
        progressbar(ijk/num_roots);
%
%  Eigenvectors & Participation Factors
%
        Dc(ijk)=0;
        Cc(ijk)=0;
        if(root(ijk)>0.01)
            L=sqrt(root(ijk));
            Dc(ijk)=besselj(n,L)/besseli(n,L);
            A=Dc(ijk);
            beta=L/radius;
            x=beta*radius;
%
            Q=0;
            Q = quad(@circular_ss_bessel_function,0,x,[],[],n,A);
            Cc(ijk)=L/sqrt(Q);
%
            Q=0;
            if(n==0)
                Q = quad(@circular_ss_bessel_function_PF,0,x,[],[],n,A);
            end
            PF(ijk)=2*Cc(ijk)*Q/L^2;
%
        end
%
        if(k==nm1)
            break;
        end
        k=k+1;
        ijk=ijk+1;
%        
    end 
end
%
progressbar(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% plot(x,y);
% grid on;
% axis([0,max(x),-10,10]);
disp(' ');
root';
%
clear RM;
clear n;
clear k;
%
sz=size(Cc);
if(sz(2)>sz(1))
    Cc=Cc';
end
sz=size(Dc);
if(sz(2)>sz(1))
    Dc=Dc';
end
sz=size(PF);
if(sz(2)>sz(1))
    PF=PF';
end
%
RM=[root' rootn' rootk' Cc Dc PF];
www=RM;
%
RM = sortrows(RM,1);
%
sz=size(root);
%
for i=max(sz):-1:1
    if(RM(i,1)<=0.001)
       RM(i,:)=[];
    end
end    
%
clear Cc;
clear Dc;
%
root=RM(:,1);
n=RM(:,2);
k=RM(:,3);
Cc=RM(:,4);
Dc=RM(:,5);
PF=RM(:,6);
%
disp(' ');
disp('  Natural ');
disp(' Frequency    ');
disp('   (Hz)    n  k       C        D        root      PF    EMM');
% 
DTerm=sqrt(D/(rho*h))/(tpi*radius^2);
clear f;
clear nn;
clear kk;
j=1;
%
sz=size(root);
NT=max(sz);
if(NT>14)
    NT=14;
end
%
fn=zeros(NT,1);
f=zeros(NT,1);
%
for i=1:NT
        fn(i)=root(i)*DTerm;
        EMM=(PF(i))^2;
        out1=sprintf('%8.2f   %d  %d   %8.5f %9.6f %8.4f %8.4f %8.4f',fn(i),n(i),k(i),Cc(i),Dc(i),sqrt(root(i)),PF(i),EMM);
        disp(out1);
        f(j)=fn(i);
        nn(j)=n(i);
        kk(j)=k(i);
        j=j+1;
end
%
nodal_diam=n;
%
disp(' ')
disp('   n = nodal diameters ');
disp('   k = nodal circles ');
disp('  PF = participation factor ');
disp(' EMM = effective modal mass '); 
%
%*********************************************************
%
fig_num=1;
%
disp(' ');
disp(' plot mode shapes? ');
disp(' 1=yes 2=no ');
ims=input(' ');
if(ims==1)
    clear n;
    mm=90;
    mmp=mm+1;
    for i=1:4
        clear Z;
        clear r;
        clear theta;
        n=nn(i);
        lambda(i)=sqrt(root(i));
%
        L=lambda(i);
%
        for ii=1:mmp
            theta(ii)=(ii-1)*2*pi/mm;
            arg=n*theta(ii);
            for(jj=1:mmp)
                r(jj)=(jj-1)/mm;  % note that the radius is non-dimensional
                kr=L*r(jj);
                Z(ii,jj)=Cc(i)*( besselj(n,kr)- Dc(i)*besseli(n,kr))*cos(arg);
            end
        end
%        
        figure(fig_num);
        fig_num=fig_num+1;
%  
        polarplot3d(Z','PlotType','surfn','meshscale',0.5,'PolarGrid',[20 72]);
%
        out1=sprintf('Mode %d  fn=%8.4g Hz  n=%d  k=%d',i,f(i),nn(i),kk(i));
        title(out1);
        h = get(gca, 'title');
        set(h, 'FontName', 'Arial','FontSize',12); 
%
    end
end
%
idamp=0;
%
disp(' ');
disp(' Calculate Base Excitation Transfer function?  1=yes  2=no ');
%
icf=input(' ');
%
if(icf==1)
%
%*********************************************************
%
%  Calculate the Base Excitation FRF
%
    disp(' ');
    damp=input(' Enter uniform damping ratio: ');
    idamp=1;
%
    disp(' ');
    f1=input('      Enter starting freq (Hz): ');
    f2=input('        Enter ending freq (Hz): ');
    df=input(' Enter the frequency step (Hz): ');
    disp(' ');
%
    num=ceil((f2-f1)/df);
%
    j=sqrt(-1);
%
    Hrd=zeros(num,1);
    Hra=zeros(num,1);
%
    freq=zeros(num,1);
%
    progressbar % Create figure and set starting time
%
    for i=1:num
        progressbar(i/num);
        freq(i)=f1+(i-1)*df;
        omega=tpi*freq(i);
%       
        for jk=1:NT
            if(abs(PF(jk))>0)
                omegan=tpi*fn(jk);
                den=(omegan^2-omega^2)+j*(2*damp*omegan*omega);
%
                kr=0; 
                arg=0.;
                nn=nodal_diam(jk);
%        
                ZZ=Cc(jk)*( besselj(nn,kr)- Dc(jk)*besseli(nn,kr))*cos(arg);  
%
                Hrd(i)=-Hrd(i)+ZZ*PF(jk)/den;
                Hra(i)= Hra(i)+ZZ*PF(jk)*omega^2/den;
            end
%                             
        end
    end
    progressbar(1);
%
    sz=size(freq);
    if(sz(2)>sz(1))
        freq=freq';
    end
%
    Haa=1+Hra;
%
    Hrd_phase=(180/pi)*atan2(imag(Hrd),real(Hrd)); 
    Haa_phase=(180/pi)*atan2(imag(Haa),real(Haa)); 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    figure(fig_num);
    fig_num=fig_num+1;
%
    if(iu==1)
        Hrd=Hrd*386.;
    else
        Hrd=Hrd*9.81;
    end
%
    Haa_mag=abs(Haa);
    Hrd_mag=abs(Hrd);
%
    max_rd=max(Hrd_mag);
    index_rd=find(Hrd_mag==max_rd);
%
    disp(' ');
    disp('  Maximum Transfer Magnitude Values ');
%
    if(iu==1)
        out1=sprintf('\n Relative Displacement = %8.4g in/G at %7.4g Hz',max_rd,freq(index_rd));
    else
        out1=sprintf('\n Relative Displacement = %8.4g m/G at %7.4g Hz',max_rd,freq(index_rd));    
    end
%
    disp(out1);
%
    max_aa=max(Haa_mag);
    index_aa=find(Haa_mag==max_aa);
%
    out1=sprintf(' Absolute Acceleration = %8.4g G/G at %7.4g Hz',max_aa,freq(index_aa));
    disp(out1);
%
    subplot(3,1,1);
    plot(freq,Hrd_phase);
    grid on;
    ylabel('Phase (deg)');
    h = get(gca, 'ylabel');set(h, 'FontName', 'Arial','FontSize',12); 
    title('Relative Displacement/Base Acceleration Transfer Function at Center of Plate ');
    axis([f1,f2,-180,180]);
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','lin','ytick',[-180,-90,0,90,180]);
    h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',12); 
    
%
    subplot(3,1,[2 3]);
    plot(freq,Hrd_mag);
    grid on;
    xlabel('Freq(Hz)'); 
%
    if(iu==1)
        ylabel('Rel Disp/Accel (in/G) ');
    else
        ylabel('Rel Disp/Accel (m/G) ');   
    end
    h = get(gca, 'xlabel');set(h, 'FontName', 'Arial','FontSize',12); 
    h = get(gca, 'ylabel');set(h, 'FontName', 'Arial','FontSize',12);     
%
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
%
    ymax= 10^(ceil(log10(max(abs(Hrd)))+0.1));
    ymin= 10^(floor(log10(min(abs(Hrd)))-0.1));
%
    LL=ymax/10000;
%
    if(ymin<LL)
        ymin=LL;
    end
%
    axis([f1,f2,ymin,ymax]);
%
    figure(fig_num);
    fig_num=fig_num+1;
%
    subplot(3,1,1);
    plot(freq,Haa_phase);
    grid on;
    title('Response Accel/Base Accel Transfer Function at Center of Plate ');
    ylabel('Phase (deg)');
    h = get(gca, 'ylabel');set(h, 'FontName', 'Arial','FontSize',12); 
    axis([f1,f2,-180,180]);
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','lin','ytick',[-180,-90,0,90,180]);
    h = get(gca, 'title');
    set(h, 'FontName', 'Arial','FontSize',12); 
%
    subplot(3,1,[2 3]);
    plot(freq,Haa_mag);
    grid on;
    xlabel('Freq(Hz)');
    ylabel('Accel/Accel (G/G) ');
    h = get(gca, 'xlabel');set(h, 'FontName', 'Arial','FontSize',12); 
    h = get(gca, 'ylabel');set(h, 'FontName', 'Arial','FontSize',12); 
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
    axis([f1,f2,min(Ylim),max(Ylim)]);
%
    accel_trans=[freq Haa_mag];
    rd_trans   =[freq Hrd_mag];
%
    accel_power_trans=[freq Haa_mag.*Haa_mag];
    rd_power_trans   =[freq Hrd_mag.*Hrd_mag];
%
    disp(' ');
    disp(' Apply sine base input?  1=yes  2=no ');
    isb=input(' ');
    if(isb==1)
%
        circular_SS_base_sine(iu,fn,damp,NT,nodal_diam,PF,Cc,Dc);
%
    end
%
end
%
disp(' ');
disp(' Apply half-sine base input?  1=yes  2=no ');
isb=input(' ');
if(isb==1)
%
   if(idamp==0)
        disp(' ');
        damp=input(' Enter uniform damping ratio: ');
        idamp=1;
   end
%
   [fig_num]=circular_SS_base_half_sine(iu,fn,damp,NT,nodal_diam,PF,Cc,Dc,fig_num);
%
end    
%