%
disp(' ');
disp(' plate_mobility.m  ver 1.0  May 18, 2011 ');
disp(' by Tom Irvine ');
disp(' ');
disp(' This program calculates the mobility of a');
disp(' simply-supported plate excited by a point force. ');
%
clear f;
clear mobility;
clear accelerance;
clear compliance;
clear omegamn;
clear omn;
clear om2;
clear s1;
clear s2;
clear s1s2;
%
disp(' ');
disp(' Enter unit type: 1=English 2=metric');
iu=input(' ');
disp(' ');
ib=1;
%
    disp(' ');
    imat=input(' Enter material: \n 1=aluminum  2=steel  3=G10  4=other ');
%
    if(imat==1) % Al
        if(iu==1)
            rho=0.1;
            E=1.0e+07;
        else
            rho=2768;
            E=6.891e+010;
        end
        mu=0.3;
    end
    if(imat==2) % steel
        if(iu==1)
            rho=0.29;
            E=3.0e+007;
        else
            rho=7887;
            E=2.067e+011;
        end
        mu=0.3;
    end
    if(imat==3) % G10
        if(iu==1)
            rho=0.065;
            E=2.7e+006;
        else
            rho=1799;
            E=1.861e+010;        
        end
        mu=0.12;
    end
    if(imat ~=1 & imat ~=2 & imat ~=3)
        if(iu==1)
            disp(' ');   
            E = input(' Enter elastic modulus (lbf/in^2) ');
            disp(' ');
            rho = input(' Enter mass density (lbm/in^3) '); 
        else
            disp(' ');   
            E = input(' Enter elastic modulus (N/m^2) ');
            disp(' ');
            rho = input(' Enter mass density (kg/m^3) ');        
        end
        disp(' ');
        mu=input(' Enter Poisson ratio ');    
    end
%
disp(' ');
loss=input(' Enter loss factor ');
damp=loss/2;
%
if(iu==1)
    a=input(' Enter    length(in) ');
    b=input(' Enter     width(in)  ');
    h=input(' Enter thickness(in)  ');    
else
    a=input(' Enter    length(m) ');
    b=input(' Enter     width(m)  ');
    h=input(' Enter thickness(m)  ');   
end
%
disp(' ');
%
if(iu==1)
    x=input(' Enter response coordinate x (in) '); 
    y=input(' Enter response coordinate y (in) ');     
else
    x=input(' Enter response coordinate x (m) '); 
    y=input(' Enter response coordinate y (m) ');        
end
%
disp(' ');
%
if(iu==1)
    xo=input(' Enter input coordinate xo (in) '); 
    yo=input(' Enter input coordinate yo (in) ');     
else
    xo=input(' Enter input coordinate xo (m) '); 
    yo=input(' Enter input coordinate yo (m) ');        
end
%
disp(' ')
num=input(' Enter number of modes ');
%
disp(' ')
maxf=input(' Enter maximum frequency (Hz) ');
%
smass=(rho*(a*b*h));
disp(' ');
if(iu==1)
    out1=sprintf(' Structural mass = %8.4g lbm',smass);
else
    out1=sprintf(' Structural mass = %8.4g kg',smass);    
end
disp(out1);
%
D=E*h^3/(12*(1.-mu^2));
%
disp(' ');
if(iu==1)
    out1=sprintf(' Plate Stiffness Factor D = %8.4g (lbf in) ',D);
else
    out1=sprintf(' Plate Stiffness Factor D = %8.4g (N m) ',D);    
end
disp(out1);
%
sq_mass=sqrt(smass);
%
if(iu==1)
    rho=rho/386;
end
DD=sqrt(D/(rho*h));
%
i=1;
%
omegamn=zeros(num,num);
omn=zeros(num,num);
%
s1=zeros(num,num);
s2=zeros(num,num);
s1s2=zeros(num,num);
%
for(m=1:num)
    for(n=1:num)
%        
        omegamn(m,n)=DD*( (m*pi/a)^2 + (n*pi/b)^2 );
            omn(m,n)=omegamn(m,n);
           omn2(m,n)=(omn(m,n))^2;
%
          s1(m,n)= sin(m*pi*x/a)*sin(n*pi*y/b);
          s2(m,n)=sin(m*pi*xo/a)*sin(n*pi*yo/b);    
        s1s2(m,n)=s1(m,n)*s2(m,n);
%
        faux(i)=omegamn(m,n)/(2*pi);
        i=i+1;
    end
end
faux=sort(faux);
%
i=sqrt(-1);
df=faux(1)/200;
%
nf=round(maxf/df);
mobility=zeros(nf,1);
compliance=zeros(nf,1);
%
den=rho*h*a*b/4;
%
progressbar
for(nk=1:nf)
%    
    progressbar(nk/nf) % Update figure    
%    
    f(nk)=nk*df;
    omega=2*pi*f(nk);
    om2=omega^2;
%
    sss=0;
    for(m=1:num)
        for(n=1:num)
           omn_diff=(omn2(m,n)-om2) + i*loss*omega*omn(m,n);
           dio=1/omn_diff;
           sss=sss+dio*s1s2(m,n);      
        end
    end  
    mobility(nk)=(i*omega/den)*sss;
    accelerance(nk)=mobility(nk)*omega;
    compliance(nk)=mobility(nk)/omega;
%    
end
%
figure(1);
plot(f,abs(mobility));
%
out1=sprintf('Point Mobility');
title(out1);
xlabel('Frequency (Hz)');
if(iu==1)
    ylabel('Velocity/Force (in/sec)/lbf');
else
    ylabel('Velocity/Force (m/sec)/N');    
end    
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','lin','YScale','log');
grid on;
%
figure(2);
plot(f,abs(accelerance));
%
out1=sprintf('Point Accelerance');
title(out1);
xlabel('Frequency (Hz)');
if(iu==1)
    ylabel('Acceleration/Force (in/sec^2)/lbf');
else
    ylabel('Acceleration/Force (m/sec^2)/N');    
end    
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','lin','YScale','log');
grid on;
%
figure(3);
plot(f,abs(compliance));
%
out1=sprintf('Point Compliance');
title(out1);
xlabel('Frequency (Hz)');
if(iu==1)
    ylabel('Displacement/Force (in/lbf)');
else
    ylabel('Displacement/Force (m/N)');    
end    
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','lin','YScale','log');
grid on;
%
CL=sqrt(E/rho);
if(iu==1)
    c=13500;
else
    c=343;
end
cr=c^2/(1.8*CL*h);
disp(' ');
out1=sprintf(' Critical Frequency in air = %6.4g Hz ',cr);
disp(out1);
%
Z=2.3*CL*rho*h^2;
tc=1/Z;
disp(' ');
if(iu==1)
    out1=sprintf(' Theoretical Compliance = %6.4g (in/sec)/lbf ',tc);
else
    out1=sprintf(' Theoretical Compliance = %6.4g (m/sec)/N  ',tc);  
end
disp(out1);
%