%
%  ss_ss_ss_ss_plate_uniform_altb.m  ver 1.0   by Tom Irvine
%
function[disp_transfer,von_Mises_stress_transfer,f,...
        HM_stress_xx,HM_stress_yy,HM_stress_xy,...
        H_stress_xx,H_stress_yy,H_stress_xy,fbig]=...
ss_ss_ss_ss_plate_uniform_altb(E,D,rho,h,a,b,mu,mass_per_area,damp,fmin,fmax,iu,x,y)
%
fig_num=2;

nmodes=6;
%
mass=mass_per_area*a*b;
%
DD=sqrt(D/(rho*h));
%
omegamn=zeros(nmodes,nmodes);
faux=zeros(nmodes*nmodes,1);

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(mass);
%
AAA=(2*sqrt(mass)/pi^2);

part=zeros(nmodes,nmodes);
fbig=zeros(nmodes*nmodes,6);

iv=1;

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);
mt=iv-1;
%
% x=a/2;
% y=b/2;

disp(' ');
disp('    fn(Hz)   m   n        PF       EMM     EMM ratio   Eig at (x,y)');
%
fn=zeros(mt,1);
for i=1:mt  
    m=fbig(i,2);
    n=fbig(i,3);
    SS=sin(m*pi*x/a)*sin(n*pi*y/b);
    eig=Amn*SS;
    fbig(i,6)=eig;
    out1=sprintf(' %9.5g \t %d\t %d\t %8.4g\t %8.4g\t %8.4g\t %8.4g  ',fbig(i,1),fbig(i,2),fbig(i,3),fbig(i,4),fbig(i,5),fbig(i,5)/mass,eig);
    disp(out1);
    fn(i)=fbig(i,1);
end
disp(' ');
if(x>a)
   disp(' Warning: x reset to total length');
   x=a;
end
%
disp(' ');
if(y>b)
   disp(' Warning: y reset to total length');
   y=b;
end
%
nf=20000;
%
clear f;
f(1)=fmin;
for k=2:nf
    f(k)=f(k-1)*2^(1/48);
    if(f(k)>fmax)
        break;
    end    
end
nf=max(size(f));
%
clear omega;
clear HM;
clear HM_stress_xx;
clear HM_stress_yy;
clear HM_stress_xy;
clear HM_stress_vM;
%
omega=2*pi*f;
%%
j=sqrt(-1); 
%
H=zeros(nf,1);
H_stress_xx=zeros(nf,1);  
H_stress_yy=zeros(nf,1);
H_stress_xy=zeros(nf,1);
HM_stress_vM=zeros(nf,1);

for k=1:nf
%    
    for m=1:nmodes
        for n=1:nmodes
%
            CC=cos(m*pi*x/a)*cos(n*pi*y/b);
            SS=sin(m*pi*x/a)*sin(n*pi*y/b);
%            
            G=(cos(m*pi)-1)*(cos(n*pi)-1);
%            
            num=G*SS;
%            
            den=(omegamn(m,n)^2-omega(k)^2) + j*2*damp*omega(k)*omegamn(m,n);
            den=den*(m*n);
%
            H(k)=H(k)+num/den;
%
            dxx=-G*(m*pi/a)^2*SS;
            dyy=-G*(n*pi/b)^2*SS;
            dxy= G*(m*pi/a)*(n*pi/b)*CC;
%
            H_stress_xx(k)=H_stress_xx(k)+ (dxx + mu*dyy)/den;
            H_stress_yy(k)=H_stress_yy(k)+ (mu*dxx + dyy)/den;
            H_stress_xy(k)=H_stress_xy(k)+  dxy/den;
%
        end
    end
end
%
P=4/(rho*h*pi^2);
%
H=P*H;
H_stress_xx=P*H_stress_xx;
H_stress_yy=P*H_stress_yy;
H_stress_xy=P*H_stress_xy;
%
zh=h/2;
AA=-E*zh/(1-mu^2);
BB=-E*zh/(1+mu);
%
H_stress_xx=H_stress_xx*AA;
H_stress_yy=H_stress_yy*AA;
H_stress_xy=H_stress_xy*BB;
%
A=[1 -0.5 0; -0.5 1 0; 0 0 3];

for k=1:nf
    sxx=H_stress_xx(k);   
    syy=H_stress_yy(k);   
    sxy=H_stress_xy(k);  

    sss=[sxx; syy; sxy];

    vvv=(conj(sss))'*A*sss;

    HM_stress_vM(k)=sqrt(abs(vvv));
end
%
%%%

HM=abs(H);
HM_stress_xx=abs(H_stress_xx);
HM_stress_yy=abs(H_stress_yy);
HM_stress_xy=abs(H_stress_xy);

[fig_num]=acoustic_plate_plots_alt(fig_num,x,y,iu,f,H,H_stress_xx,H_stress_yy,H_stress_xy,HM_stress_vM);
%%%

%
f=fix_size(f);
HM=fix_size(HM);

LL=min([length(f) length(HM)]);

xc=[f(1:LL) HM(1:LL)];
disp_transfer=xc;
clear C;
clear I;
[C,I]=max(xc);
xmax=xc(I(2),2);
fmax=xc(I(2),1);
%
if(iu==1)
    out1=sprintf('Displacement: max = %8.4g in/psi at %8.4g Hz  ',xmax,fmax);
else
    out1=sprintf('Displacement: max = %8.4g m/Pa at %8.4g Hz  ',xmax,fmax);    
end
disp(out1);
%
%
%%%%%%
%
HM_stress_vM=fix_size(HM_stress_vM);

xc=[f(1:LL) HM_stress_vM(1:LL)];
von_Mises_stress_transfer=xc;

clear C;
clear I;
[C,I]=max(xc);
xmax=xc(I(2),2);
fmax=xc(I(2),1);
%
disp(' ');
if(iu==1)
    out1=sprintf('von Mises stress: max = %8.4g (psi/psi) at %8.4g Hz  ',xmax,fmax);
else
    out1=sprintf('von Mises stress: max = %8.4g (Pa/Pa) at %8.4g Hz  ',xmax,fmax);    
end
disp(out1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear xx;
clear yy;
clear zzr;
%
num=41; 
%
clear x;
clear y;
%
dx=a/(num-1);
dy=b/(num-1);
%
x=zeros(num,1);
y=zeros(num,1);

for i=1:num
   ii=i-1;
   x(i)=ii*dx;
   y(i)=ii*dy;
end
%
%  Calculate mode shapes and mass-normalized
%
clear norm;
%
zzr=zeros(num,num);
%
for i=1:num 
%          
        for j=1:num 
%
            zzr(i,j)=sin(pi*x(i)/a)*sin(pi*y(j)/b);            
%
        end
end

%
norm=sqrt(mass_per_area*a*b);
%
zzr=2*zzr/norm;
%
zmax=max(max(zzr));
zmin=min(min(zzr));
%
clear aaa;
clear abc;
abc=[x y];
aaa=max(abc);
xmin=0;
ymin=0;
xmax=aaa;
ymax=aaa;
figure(1); 
surf(x,y,zzr);
zmax=2*zmax;
% axis([xmin,xmax,ymin,ymax,zmin,zmax]);
zlim([zmin zmax]);
% ylim([0 b]);
% xlim([0 a]);
out1=sprintf('Mode 1  fn=%8.4g Hz',fn(1));
title(out1)
if(iu==1)
    xlabel('X (in)');
    ylabel('Y (in)');    
else
    xlabel('X (m)');
    ylabel('Y (m)');
end    