disp(' ');
disp(' ss_plate_base.m  ver 1.1  August 19, 2009');
disp(' by Tom Irvine  Email: tomirvine@aol.com ');
disp(' ');
disp(' Normal Modes & Optional Base Excitation ');
disp(' ');
%
clear root;
clear f;
clear H;
clear HA;
clear HA2;
clear beta;
clear fn;
clear fmn;
clear faux;
clear A;
clear M;
clear D;
clear x;
clear y;
clear z;
clear YY;
clear omegamn;
clear fbig;
clear part;
%
nmodes=5;
mt=nmodes^2;
%
disp(' ');
%
length=input(' Enter the length (inch) ');
width=input(' Enter the width (inch) ');
thick=input(' Enter the thickness (inch) ');
A=width*thick;
%
a=length;
b=width;
h=thick;
%
disp(' ');
imat=input(' Enter material: \n 1=aluminum  2=steel  3=G10  4=other ');
%
if(imat==1) % Al
    rho=0.1;
	E=1.0e+07;
    mu=0.3;
end
if(imat==2) % steel
    rho=0.29;
    E=3.0e+007;
    mu=0.3;
end
if(imat==3) % G10
    rho=0.065 ;
    E=2.7e+006;
    mu=0.12;
end
if(imat ~=1 & imat ~=2 & imat ~=3)
    disp(' ');   
    E = input(' Enter the elastic modulus (lbf/in^2) ');
    disp(' ');
    rho = input(' Enter the mass density (lbm/in^3) '); 
    disp(' ');
    mu=input(' Enter Poisson ratio ');    
 end
%
smass=(rho*(a*b*h));
disp(' ');
out1=sprintf(' Structural mass = %8.4g lbm',smass);
disp(out1);
disp(' ');
disp(' Add non structural mass ?');
disp(' 1=yes  2=no ')
imass=input(' ');
%
NSM=0;
if(imass==1)
  disp(' ');
  disp(' Enter added mass(lbm) ');
  NSM=input(' ');
end
%
mass=NSM+smass;
M=mass;
%
disp(' ');
out1=sprintf(' Total mass = %8.4g lbm ',M);
disp(out1);
%
mass=mass/386;
M=M/386;
rho=rho/386;   % mass per volume
%
D=E*h^3/(12*(1.-mu^2));
%
disp(' ');
out1=sprintf(' Plate Stiffness Factor D = %8.4g (lbf in) ',D);
disp(out1);
%
sq_mass=sqrt(mass);
%
DD=sqrt(D/(rho*h));
%
i=1;
for(m=1:nmodes)
    for(n=1:nmodes)
        omegamn(m,n)=DD*( (m*pi/a)^2 + (n*pi/b)^2 );
        faux(i)=omegamn(m,n)/(2*pi);
        i=i+1;
    end
end
%
fmn=omegamn/(2*pi);
sort(faux);
%
Amn=2/sqrt(M);
%
AAA=(2*sqrt(M)/pi^2);
iv=1;
part=zeros(nmodes,nmodes);
for(i=1:nmodes)
    for(j=1:nmodes)
        part(i,j)=(cos(i*pi)-1)*(cos(j*pi)-1);
        part(i,j)=AAA*part(i,j)/(i*j);
        fbig(iv,1)=faux(iv);
        fbig(iv,2)=i;
        fbig(iv,3)=j;
        fbig(iv,4)=part(i,j);
        fbig(iv,5)=(part(i,j))^2;
        iv=iv+1;
    end
end
fbig=sortrows(fbig,1);
%
disp(' ');
disp('    fn(Hz)   m   n        PF    EMM ratio');
for(i=1:mt)  
    out1=sprintf(' %9.5g \t %d\t %d\t %8.4g\t %8.4g  ',fbig(i,1),fbig(i,2),fbig(i,3),fbig(i,4),fbig(i,5)/M);
    disp(out1);
end
disp(' ');
%
dx=a/200;
dy=b/200;
%
pa=pi/a;
pb=pi/b;
%
clear z;
clear zmn;
%
x=zeros(1,201);
y=zeros(1,201);
z=zeros(201,201);
zmn=zeros(mt,201,201);
%
for(i=1:201)
    x(i)=pa*(i-1)*dx;   
    y(i)=pb*(i-1)*dy;  
end
%
for(iv=1:mt)
    for(i=1:201)
        term1=sin(fbig(iv,2)*x(i));
        for(j=1:201)    
            z(i,j)=term1*sin(fbig(iv,3)*y(j));
            zmn(iv,i,j)=z(i,j);
        end        
    end
    if(iv==1)
        figure(1);
        zzr=z';
        surf(x,y,zzr);
        out1=sprintf('First Mode  fn=%8.4g Hz',fbig(1,1));
        title(out1);
    end
    if(iv==2)
        figure(2); 
        zzr=z';
        surf(x,y,zzr);
        out1=sprintf('Second Mode  fn=%8.4g Hz',fbig(2,1)); 
        title(out1); 
    end  
end
%
z=z*Amn;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Calculate Frequency Response Function ');
disp(' 1=yes  2=no ');
%
ia=input(' ');
%
if(ia==1)
%
    clear f;
    clear x;
    clear y;
    clear damp;
    clear Z;
    clear H;
    clear HA;
    clear om;
%    
    disp(' ');    
    out1=sprintf(' Enter uniformal modal damping ratio ');
    disp(out1);
    damp=input(' ');
%    
    disp(' ');
    disp(' Enter distance x ');
    x=input(' ');
    if(x>a)
        disp(' Warning: x reset to total length');
        x=a;
    end
%
    disp(' ');
    disp(' Enter distance y ');
    y=input(' ');
    if(y>b)
        disp(' Warning: y reset to total length');
        y=b;
    end
%
   disp(' ');
   disp(' Enter maximum base excitation frequency Hz ');
   MAXF=input(' ');
%
   nf=20000;
   f(1)=1;
   for(k=2:nf)
        f(k)=f(k-1)*2^(1/48);
        if(f(k)>MAXF)
            break;
        end    
   end
   nf=max(size(f));
%
    pax=x*pi/a;
    pby=y*pi/b;
%
   Amn=2/sqrt(M);
%
    sq_mass=sqrt(M);
%
    tpi=2*pi;
    j=sqrt(-1);   
%
    clear H;
    clear HA;
%
    for(k=1:nf)
        H(k)=0;
        HA(k)=0;
        om=tpi*f(k);
 %
        for(i=1:nmodes)
            for(jk=1:nmodes)          
 %                  
                clear num;
                clear den;
                clear pY;
                clear omn;
 %               
                nmode=Amn*sin(i*pax)*sin(jk*pby);  
                pY=part(i,jk)*nmode;
                omn=omegamn(i,jk);
                num=-pY;
                den=(omn^2-om^2)+j*2*damp*om*omn;
                H(k)=H(k)+num/den;
 %
                num=om^2*pY;
                HA(k)=HA(k)+(num/den);
 %      
            end
 %
        end
    end
%
    H=386*abs(H);
    HA=abs(HA)+1;
    HA2=HA.*HA;
%
    maxH=0.;
    maxHA=0.;
%
    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       
%
    figure(3);
    plot(f,H);
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
    out1=sprintf('Rel Disp Frequency Response Function  x=%7.4g in  y=%7.4g in',x,y);
    title(out1);
    ylabel('Rel Disp (in)/ Base Accel (G) ');
    xlabel('Frequency (Hz)');
    grid on;
%
    figure(4);
    plot(f,HA);
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
    out1=sprintf('Accel Frequency Response Function  x=%7.4g in  y=%7.4g in',x,y);
    title(out1);
    ylabel('Response Accel / Base Accel ');
    xlabel('Frequency (Hz)');
    grid on; 
%    
    figure(5);
    plot(f,HA2);
    set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
    out1=sprintf('Power Transmissibility  x=%7.4g in  y=%7.4g in',x,y);
    title(out1);
    ylabel('Trans (G^2/G^2) ');
    xlabel('Frequency (Hz)');
    grid on;
    disp(' ');
%
    disp(' ');
    out1=sprintf('  max Rel Disp FRF = %8.3e (in/G) at %8.4g Hz ',maxH,maxF);
    disp(out1);
    disp(' ');
    out1=sprintf('  max Accel FRF    = %8.4g (G/G)     at %8.4g Hz ',maxHA,maxFA);
    disp(out1);
    out1=sprintf('  max Power Trans  = %8.4g (G^2/G^2) at %8.4g Hz ',maxHA^2,maxFA);
    disp(out1); 
%%
end