function test

% Konstanten
lx  = 2.8; 
ly  = 2; 
h   = 0.01; 
rho = 2.3e3; 
eta = 0.01; 
x1  = 1.99;
x2  = 1.99; 
y1  = 1.27; 
y2  = 1.27; 
v   = 0.2; 
E   = 2.6e10;

% m und n als Vektoren
m = 1:20; 
n = 1:20;

% Einzelterme
[Psi1mn,Psi2mn] = Psi(m,n);
Psimn = Psi1mn.*Psi2mn; 

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

YWf = zeros(5000,1); f = 1:5000;

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

loglog(f,YWf)
grid on;

    function [Psi1mn,Psi2mn] = Psi(m,n)

        yim1 = 0.5*pi*(4*m-1); yim1(m==1) = 4.73004;
        yjm2 = 0.5*pi*(4*m+1); yjm2(m==1) = 7.8532;
        yin1 = 0.5*pi*(4*n-1); yin1(n==1) = 4.73004;
        yjn2 = 0.5*pi*(4*n+1); yjn2(n==1) = 7.8532;

        mu = rem(m,2)~=0; % ungerade Elemente von m
        mg = rem(m,2)==0; % gerade Elemente von m
        nu = rem(n,2)~=0;
        ng = rem(n,2)==0;
        
        phix1m = zeros(size(m));
        phiy1n = zeros(size(n));
        
        phix1m(mu) = sqrt(2)* (cos(yim1(mu)) * (x1/lx-0.5) - ...
            sin(yim1(mu)*0.5) ./ sinh(yim1(mu)*0.5) .* cosh(yim1(mu)) * (x1/lx-0.5));
        phix1m(mg) = sqrt(2)* (cos(yjm2(mg)) * (x1/lx-0.5) + ...
            sin(yjm2(mg)*0.5) ./ sinh(yjm2(mg)*0.5) .* sinh(yjm2(mg)) * (x1/lx-0.5));

        phiy1n(nu) = sqrt(2) * (cos(yin1(nu)) * (y1/ly-0.5) - ...
            sin(yin1(nu)*0.5) ./ sinh(yin1(nu)*0.5) .* cosh(yin1(nu)) * (y1/ly-0.5));
        phiy1n(ng) = sqrt(2) * (cos(yjn2(ng)) * (y1/ly-0.5) + ...
            sin(yjn2(ng)*0.5) ./ sinh(yjn2(ng)*0.5) .* sinh(yjn2(ng)) * (y1/ly-0.5));

        Psi1mn = phix1m'*phiy1n;
        
        phix2m = zeros(size(m));
        phiy2n = zeros(size(n));
       
        phix2m(mu) = sqrt(2) * (cos(yim1(mu)) * (x2/lx-0.5) - ...
            sin(yim1(mu)*0.5) ./ sinh(yim1(mu)*0.5) .* cosh(yim1(mu)) * (x2/lx-0.5));
        phix2m(mg) = sqrt(2) * (cos(yjm2(mg)) * (x2/lx-0.5) + ...
            sin(yjm2(mg)*0.5) ./ sinh(yjm2(mg)*0.5) .* sinh(yjm2(mg)) * (x2/lx-0.5));
        
        phiy2n(nu) = sqrt(2) * (cos(yin1(nu)) * (y2/ly-0.5) - ...
            sin(yin1(nu)*0.5) ./ sinh(yin1(nu)*0.5) .* cosh(yin1(nu)) * (y2/ly-0.5));
        phiy2n(ng) = sqrt(2) * (cos(yjn2(ng)) * (y2./ly-0.5) + ...
            sin(yjn2(ng)*0.5) ./ sinh(yjn2(ng)*0.5) .* sinh(yjn2(ng)) * (y2/ly-0.5));
                
        Psi2mn = phix2m'*phiy2n;        
    end % function Psi

    function wmn = w(m,n)
        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 w

end % function test