function YWFZ

% constants
lx  = 2.8;    % x-length
ly  = 2;      % y-length
h   = 0.1;    % thickness of the plate
rho = 2.3e3;  % Density
eta = 0.03;   % loss factor
x1  = 1.99;
x2  = 1.99; 
y1  = 1.27;  
y2  = 1.27; 
v   = 0.2;    % Poisson's ratio
E   = 2.6e10; % Young's modulus

% m und n as vectors
m = 1:20; % bending mode
n = 1:20; % bending mode

% natural mode
[Psi1mn,Psi2mn] = Psi(m,n);
Psimn = Psi1mn'*Psi2mn; 

% natural frequency
w2mn = w(m,n).^2;

% mobility
f = 1:3000;
YWf = zeros(length(f),1); 

    for F=f;
    matrix_mn = Psimn ./ ((1+j*eta) * w2mn - (2*pi*F)^2);
    Summe  = sum(matrix_mn(:));
    YWf(F)  = abs(((2*pi*j*F) / (rho*lx*ly*h)) * Summe);
    end

    loglog(f,YWf)
    grid on;

    %Characteristic beam function
    
    function [Psi1mn,Psi2mn] = Psi(m,n)
      
        phix1m=zeros(length(m)); 
        phiy1n=zeros(length(n));
    
       for M = m
           
           if rem(M,2)
           ii = (M+1)/2;
           yim1 = 0.5*pi*(4*ii-1); yim1(ii==1) = 4.73004;
           phix1m = sqrt(2)*(cos(yim1)*(x1/lx-0.5)-...
           sin(yim1*0.5)./sinh(yim1*0.5).*cosh(yim1)*(x1./lx-0.5)); 
           
           else
           jj = M/2;
           yjm2 = 0.5*pi*(4*jj+1); yjm2(jj==1) = 7.8532;
           phix1m =sqrt(2)*(cos(yjm2)*(x1/lx-0.5)+...
           sin(yjm2*0.5)./sinh(yjm2*0.5).*sinh(yjm2)*(x1./lx-0.5)); 
           end 
       end
       
       for N = n
        
           if rem(N,2)
           ii = (N+1)/2;
           yin1 = 0.5*pi*(4*ii-1); yin1(ii==1) = 4.73004;
           phiy1n=sqrt(2)*(cos(yin1)*(y1/ly-0.5)-...
           sin(yin1*0.5)./sinh(yin1*0.5).*cosh(yin1)*(y1./ly-0.5));
       
           else
           jj = N/2;
           yjn2 = 0.5*pi*(4*jj+1); yjn2(jj==1) = 7.8532;
           phiy1n=sqrt(2)*(cos(yjn2)*(y1/ly-0.5)+...
           sin(yjn2*0.5)./sinh(yjn2*0.5).*sinh(yjn2)*(y1./ly-0.5));
           end
       end   

    Psi1mn = phix1m'*phiy1n;
        
        phix2m=zeros(length(m)); 
        phiy2n=zeros(length(n));
        
       for M = m
        
        if rem(M,2)
        ii = (M+1)/2; 
        yim1 = 0.5*pi*(4*ii-1); yim1(ii==1) = 4.73004;
        phix2m=sqrt(2)*(cos(yim1)*(x2/lx-0.5)-...
        sin(yim1*0.5)./sinh(yim1*0.5).*cosh(yim1)*(x2/lx-0.5)); 
        else
        jj = M/2;
        yjm2 = 0.5*pi*(4*jj+1); yjm2(jj==1) = 7.8532;
        phix2m=sqrt(2)*(cos(yjm2)*(x2/lx-0.5)+...
        sin(yjm2*0.5)./sinh(yjm2*0.5).*sinh(yjm2)*(x2/lx-0.5));
        end
       end

       for N = n
    
        if rem(N,2)
        ii = (N+1)/2;
        yin1 = 0.5*pi*(4*ii-1); yin1(ii==1) = 4.73004;
        phiy2n=sqrt(2)*(cos(yin1)*(y2/ly-0.5)-...
        sin(yin1*0.5)./sinh(yin1*0.5).*cosh(yin1)*(y2/ly-0.5));
        else
        jj = N/2;
        yjn2 = 0.5*pi*(4*jj+1); yjn2(jj==1) = 7.8532;
        phiy2n=sqrt(2)*(cos(yjn2)*(y2/ly-0.5)+...
        sin(yjn2*0.5)./sinh(yjn2*0.5).*sinh(yjn2)*(y2/ly-0.5));
        end
       end
    
                
   Psi2mn = phix2m'*phiy2n;        
    
    end % function Psi

    function wmn = w(m,n)
        
        %Boundary Conditions
        
        Gxm = m + 0.5; Gxm(m==1) = 1.506;

        Gyn = n + 0.5; Gyn(n==1) = 1.506;

        Hxm = (m+0.5).^2 .*(1-4./(pi*(2*m+1)));
        Hxm(m==1) = 1.248;

        Hyn = (n+0.5).^2 .*(1-4./(pi*(2*n+1)));
        Hyn(n==1) = 1.248;

        Jxm = (m+0.5).^2 .*(1+12./(pi*(2*m+1)));
        Jxm(m==1) = 5.017;

        Jyn = (n+0.5).^2 .*(1+12./(pi*(2*n+1)));
        Jyn(n==1) = 5.017;

        Hmn = Hxm'*Hyn;
        Jmn = Jxm'*Jyn;
        
        Gxm_matrix = ones(length(m),1)*Gxm;
        Gyn_matrix = Gyn'*ones(1,length(n));     
                
        qmn = sqrt(Gxm_matrix.^4 + Gyn_matrix.^4 .* (lx/ly)^4+2*(lx/ly)^2 ...
              * (v*Hmn + (1-v)*Jmn));

        wmn = sqrt(E*h^2/(12*rho*(1-v^2))) * (pi/lx)^2 * qmn;

    end % function wmn

end % YWFZ