clear all;

% constants
f = linspace(1.25,4000,3200);
roh = 2.3e3; % Density                           
E = 2.6e10;  % Young's modulus                         
eta = 0.01;  % loss factor                     
h = 0.10;    % thickness of the plate                        
nue = 0.3;   % Poisson's ratio
lx = 2.8;                                
ly = 2;   
x1 = 1.99;                            
y1 = 1.27;                              
x2 = x1;                              
y2 = y1;
x3 = 1;
y3 = 1;
x4 = 1;
y4 = 1;
x5 = x1;
y5 = y1;
x6 = 1;
y6 = 1;
x7 = 1;
y7 = 1;
x8 = 1;
y8 = 1;
w=(2*pi).*f;                        

% max. frequency
fmax=5000;               

%zeros of the gamma-functions
yi = [7.8532 14.13716];
yj = [4.73004 10.9956];

groessef=size(f);
zz = zeros(1,groessef(2)); 
Yf1z1zsummeges{15} = zz;
Yf1z1zsummeges{16} = zz;
Yf1z1zsummeges{17} = zz;
Yf1z1zsummeges{18} = zz;
Yf1z1zsummeges{25} = zz;
Yf1z1zsummeges{26} = zz;
Yf1z1zsummeges{27} = zz;
Yf1z1zsummeges{28} = zz;
Yf1z1zsummeges{35} = zz;
Yf1z1zsummeges{36} = zz;
Yf1z1zsummeges{37} = zz;
Yf1z1zsummeges{38} = zz;
Yf1z1zsummeges{45} = zz;
Yf1z1zsummeges{46} = zz;
Yf1z1zsummeges{47} = zz;
Yf1z1zsummeges{48} = zz;

Yf1z1zsumme{15} = zz;
Yf1z1zsumme{16} = zz;
Yf1z1zsumme{17} = zz;
Yf1z1zsumme{18} = zz;
Yf1z1zsumme{25} = zz;
Yf1z1zsumme{26} = zz;
Yf1z1zsumme{27} = zz;
Yf1z1zsumme{28} = zz;
Yf1z1zsumme{35} = zz;
Yf1z1zsumme{36} = zz;
Yf1z1zsumme{37} = zz;
Yf1z1zsumme{38} = zz;
Yf1z1zsumme{45} = zz;
Yf1z1zsumme{46} = zz;
Yf1z1zsumme{47} = zz;
Yf1z1zsumme{48} = zz;

m = 1;
n = 0;
z = 0; %Zähler 
while (1)
    z = z + 1;
    if m == 0
        Gx = 0;
        Hx = 0;
        Jx = 0;
    elseif m == 1
        Gx = 1.506;
        Hx = 1.248;
        Jx = 5.017;
    else
        Gx = m+0.5;
        Hx = (m+0.5)^2*(1-4/((2*m+1)*pi));
        Jx = (m+0.5)^2*(1+12/((2*m+1)*pi));
    end
    if n == 0
        Gy = 0;
        Hy = 0;
        Jy = 0;
    elseif n == 1
        Gy = 1.506;
        Hy = 1.248;
        Jy = 5.017;
    else
        Gy = n+0.5;
        Hy = (n+0.5)^2*(1-4/((2*n+1)*pi));
        Jy = (n+0.5)^2*(1+12/((2*n+1)*pi));
    end
    qmn = sqrt(Gx^4+Gy^4*(lx/ly)^4+2*(lx/ly)^2*(nue*Hx*Hy+(1-nue)*Jx*Jy)); 
    wmnr = sqrt(E*h^2/(12*roh*(1-nue^2)))*(pi/lx)^2*qmn; 
    wmn = wmnr*sqrt(1+j*eta); 
    fmn = wmn/(2*pi);                                     
     
       if m==0
           phix1 = 1;
           phix2 = 1;
           phix3 = 1;
           phix4 = 1;
           phix5 = 1;
           phix6 = 1;
           phix7 = 1;
           phix8 = 1;
           
       else
           if m/2 == round(m/2); %Abfrage: m gerade Zahl
                ii = m/2;
                if ii > 2;
                   yi(ii) = (4*ii+1)*pi/2;
                end                
                kn = sin(0.5*yi(ii))/sinh(0.5*yi(ii));
                phix1 = sqrt(2)*(sin(yi(ii)*(x1/lx-0.5))+kn*sinh(yi(ii)*(x1/lx-0.5))); %phi an der Stelle x1
                phix2 = sqrt(2)*(sin(yi(ii)*(x2/lx-0.5))+kn*sinh(yi(ii)*(x2/lx-0.5))); %phi an der Stelle x2
                phix3 = sqrt(2)*(sin(yi(ii)*(x3/lx-0.5))+kn*sinh(yi(ii)*(x3/lx-0.5))); %phi an der Stelle x1
                phix4 = sqrt(2)*(sin(yi(ii)*(x4/lx-0.5))+kn*sinh(yi(ii)*(x4/lx-0.5))); %phi an der Stelle x2
                phix5 = sqrt(2)*(sin(yi(ii)*(x5/lx-0.5))+kn*sinh(yi(ii)*(x5/lx-0.5))); %phi an der Stelle x1
                phix6 = sqrt(2)*(sin(yi(ii)*(x6/lx-0.5))+kn*sinh(yi(ii)*(x6/lx-0.5))); %phi an der Stelle x2
                phix7 = sqrt(2)*(sin(yi(ii)*(x7/lx-0.5))+kn*sinh(yi(ii)*(x7/lx-0.5))); %phi an der Stelle x1
                phix8 = sqrt(2)*(sin(yi(ii)*(x8/lx-0.5))+kn*sinh(yi(ii)*(x8/lx-0.5))); %phi an der Stelle x2
                
           else % m ungerade Zahl
                jj = (m+1)/2;
                if jj > 2;
                   yj(jj) = (4*jj-1)*pi/2;
                end
                kn = -sin(0.5*yj(jj))/sinh(0.5*yj(jj));
                phix1 = sqrt(2)*(cos(yj(jj)*(x1/lx-0.5))+kn*cosh(yj(jj)*(x1/lx-0.5))); 
                phix2 = sqrt(2)*(cos(yj(jj)*(x2/lx-0.5))+kn*cosh(yj(jj)*(x2/lx-0.5)));
                phix3 = sqrt(2)*(cos(yj(jj)*(x3/lx-0.5))+kn*cosh(yj(jj)*(x3/lx-0.5))); 
                phix4 = sqrt(2)*(cos(yj(jj)*(x4/lx-0.5))+kn*cosh(yj(jj)*(x4/lx-0.5)));
                phix5 = sqrt(2)*(cos(yj(jj)*(x5/lx-0.5))+kn*cosh(yj(jj)*(x5/lx-0.5))); 
                phix6 = sqrt(2)*(cos(yj(jj)*(x6/lx-0.5))+kn*cosh(yj(jj)*(x6/lx-0.5)));
                phix7 = sqrt(2)*(cos(yj(jj)*(x7/lx-0.5))+kn*cosh(yj(jj)*(x7/lx-0.5))); 
                phix8 = sqrt(2)*(cos(yj(jj)*(x8/lx-0.5))+kn*cosh(yj(jj)*(x8/lx-0.5)));
                
           end % der if Abfrage m gerade Zahl
        end
        if n == 0
            phiy1 = 1;
            phiy2 = 1;
            phiy3 = 1;
            phiy4 = 1;
            phiy5 = 1;
            phiy6 = 1;
            phiy7 = 1;
            phiy8 = 1;
            
        else
            if n/2 == round(n/2); %Abfrage: n gerade Zahl
                ii = n/2;
                if ii > 2;
                   yi(ii) = (4*ii+1)*pi/2;
                end                
                kn = sin(0.5*yi(ii))/sinh(0.5*yi(ii));
                phiy1 = sqrt(2)*(sin(yi(ii)*(y1/ly-0.5))+kn*sinh(yi(ii)*(y1/ly-0.5))); 
                phiy2 = sqrt(2)*(sin(yi(ii)*(y2/ly-0.5))+kn*sinh(yi(ii)*(y2/ly-0.5))); 
                phiy3 = sqrt(2)*(sin(yi(ii)*(y3/ly-0.5))+kn*sinh(yi(ii)*(y3/ly-0.5))); 
                phiy4 = sqrt(2)*(sin(yi(ii)*(y4/ly-0.5))+kn*sinh(yi(ii)*(y4/ly-0.5))); 
                phiy5 = sqrt(2)*(sin(yi(ii)*(y5/ly-0.5))+kn*sinh(yi(ii)*(y5/ly-0.5))); 
                phiy6 = sqrt(2)*(sin(yi(ii)*(y6/ly-0.5))+kn*sinh(yi(ii)*(y6/ly-0.5))); 
                phiy7 = sqrt(2)*(sin(yi(ii)*(y7/ly-0.5))+kn*sinh(yi(ii)*(y7/ly-0.5))); 
                phiy8 = sqrt(2)*(sin(yi(ii)*(y8/ly-0.5))+kn*sinh(yi(ii)*(y8/ly-0.5))); 
                                
            else % n ungerade Zahl
                jj = (n+1)/2;
                if jj > 2;
                   yj(jj) = (4*jj-1)*pi/2;
                end
                kn = -sin(0.5*yj(jj))/sinh(0.5*yj(jj));
                phiy1 = sqrt(2)*(cos(yj(jj)*(y1/ly-0.5))+kn*cosh(yj(jj)*(y1/ly-0.5))); 
                phiy2 = sqrt(2)*(cos(yj(jj)*(y2/ly-0.5))+kn*cosh(yj(jj)*(y2/ly-0.5))); 
                phiy3 = sqrt(2)*(cos(yj(jj)*(y3/ly-0.5))+kn*cosh(yj(jj)*(y3/ly-0.5))); 
                phiy4 = sqrt(2)*(cos(yj(jj)*(y4/ly-0.5))+kn*cosh(yj(jj)*(y4/ly-0.5))); 
                phiy5 = sqrt(2)*(cos(yj(jj)*(y5/ly-0.5))+kn*cosh(yj(jj)*(y5/ly-0.5))); 
                phiy6 = sqrt(2)*(cos(yj(jj)*(y6/ly-0.5))+kn*cosh(yj(jj)*(y6/ly-0.5)));                
                phiy7 = sqrt(2)*(cos(yj(jj)*(y7/ly-0.5))+kn*cosh(yj(jj)*(y7/ly-0.5)));                
                phiy8 = sqrt(2)*(cos(yj(jj)*(y8/ly-0.5))+kn*cosh(yj(jj)*(y8/ly-0.5))); 
                
                
             end % der if Abfrage n gerade Zahl
        end

       psi1 = phix1*phiy1; % psi an der Stelle 1
       psi2 = phix2*phiy2; % psi an der Stelle 2
       psi3 = phix3*phiy3; % psi an der Stelle 1
       psi4 = phix4*phiy4; % psi an der Stelle 2
       psi5 = phix5*phiy5; % psi an der Stelle 1
       psi6 = phix6*phiy6; % psi an der Stelle 2
       psi7 = phix7*phiy7; % psi an der Stelle 1
       psi8 = phix8*phiy8; % psi an der Stelle 2
                                
       rr = 1./(roh.*h.*lx.*ly.*(wmn.^2-w.^2));  
       
       Yf1z1zsumme{15}(z,:) = psi1.*psi5.*rr;
       Yf1z1zsumme{16}(z,:) = psi1.*psi6.*rr;
       Yf1z1zsumme{17}(z,:) = psi1.*psi7.*rr;
       Yf1z1zsumme{18}(z,:) = psi1.*psi8.*rr;
       Yf1z1zsumme{25}(z,:) = psi2.*psi5.*rr;
       Yf1z1zsumme{26}(z,:) = psi2.*psi6.*rr;
       Yf1z1zsumme{27}(z,:) = psi2.*psi7.*rr;
       Yf1z1zsumme{28}(z,:) = psi2.*psi8.*rr;
       Yf1z1zsumme{35}(z,:) = psi3.*psi5.*rr;
       Yf1z1zsumme{36}(z,:) = psi3.*psi6.*rr;
       Yf1z1zsumme{37}(z,:) = psi3.*psi7.*rr;
       Yf1z1zsumme{38}(z,:) = psi3.*psi8.*rr;
       Yf1z1zsumme{45}(z,:) = psi4.*psi5.*rr;
       Yf1z1zsumme{46}(z,:) = psi4.*psi6.*rr;
       Yf1z1zsumme{47}(z,:) = psi4.*psi7.*rr;
       Yf1z1zsumme{48}(z,:) = psi4.*psi8.*rr;
      
       Yf1z1zsummeges{15}(1,:) = Yf1z1zsumme{15}(z,:)+Yf1z1zsummeges{15}(1,:);
       Yf1z1zsummeges{16}(1,:) = Yf1z1zsumme{16}(z,:)+Yf1z1zsummeges{16}(1,:);
       Yf1z1zsummeges{17}(1,:) = Yf1z1zsumme{17}(z,:)+Yf1z1zsummeges{17}(1,:);
       Yf1z1zsummeges{18}(1,:) = Yf1z1zsumme{18}(z,:)+Yf1z1zsummeges{18}(1,:);
       Yf1z1zsummeges{25}(1,:) = Yf1z1zsumme{25}(z,:)+Yf1z1zsummeges{25}(1,:);
       Yf1z1zsummeges{26}(1,:) = Yf1z1zsumme{26}(z,:)+Yf1z1zsummeges{26}(1,:);
       Yf1z1zsummeges{27}(1,:) = Yf1z1zsumme{27}(z,:)+Yf1z1zsummeges{27}(1,:);
       Yf1z1zsummeges{28}(1,:) = Yf1z1zsumme{28}(z,:)+Yf1z1zsummeges{28}(1,:);
       Yf1z1zsummeges{35}(1,:) = Yf1z1zsumme{35}(z,:)+Yf1z1zsummeges{35}(1,:);
       Yf1z1zsummeges{36}(1,:) = Yf1z1zsumme{36}(z,:)+Yf1z1zsummeges{36}(1,:);
       Yf1z1zsummeges{37}(1,:) = Yf1z1zsumme{37}(z,:)+Yf1z1zsummeges{37}(1,:);
       Yf1z1zsummeges{38}(1,:) = Yf1z1zsumme{38}(z,:)+Yf1z1zsummeges{38}(1,:);
       Yf1z1zsummeges{45}(1,:) = Yf1z1zsumme{45}(z,:)+Yf1z1zsummeges{45}(1,:);
       Yf1z1zsummeges{46}(1,:) = Yf1z1zsumme{46}(z,:)+Yf1z1zsummeges{46}(1,:);
       Yf1z1zsummeges{47}(1,:) = Yf1z1zsumme{47}(z,:)+Yf1z1zsummeges{47}(1,:);
       Yf1z1zsummeges{48}(1,:) = Yf1z1zsumme{48}(z,:)+Yf1z1zsummeges{48}(1,:);            
       
    if (fmn < fmax)
        m = m + 1;
    elseif (m > 0)
        n = n + 1;
        m = 0;
        z = z - 1;
    else
        break;
    end
end

xx = j.*w;

YWFz{15} = abs(Yf1z1zsummeges{15}(1,:).*xx);  % T = theta!
YWFz{16} = abs(Yf1z1zsummeges{16}(1,:).*xx);
YWFz{17} = abs(Yf1z1zsummeges{17}(1,:).*xx);
YWFz{18} = abs(Yf1z1zsummeges{18}(1,:).*xx);
YWFz{25} = abs(Yf1z1zsummeges{25}(1,:).*xx);
YWFz{26} = abs(Yf1z1zsummeges{26}(1,:).*xx);
YWFz{27} = abs(Yf1z1zsummeges{27}(1,:).*xx);
YWFz{28} = abs(Yf1z1zsummeges{28}(1,:).*xx);
YWFz{35} = abs(Yf1z1zsummeges{35}(1,:).*xx);
YWFz{36} = abs(Yf1z1zsummeges{36}(1,:).*xx);
YWFz{37} = abs(Yf1z1zsummeges{37}(1,:).*xx);
YWFz{38} = abs(Yf1z1zsummeges{38}(1,:).*xx);
YWFz{45} = abs(Yf1z1zsummeges{45}(1,:).*xx);
YWFz{46} = abs(Yf1z1zsummeges{46}(1,:).*xx);
YWFz{47} = abs(Yf1z1zsummeges{47}(1,:).*xx);
YWFz{48} = abs(Yf1z1zsummeges{48}(1,:).*xx);

Mob_Matrix1 = [YWFz{15} YWFz{16} YWFz{17} YWFz{18}
               YWFz{25} YWFz{26} YWFz{27} YWFz{28}
               YWFz{35} YWFz{36} YWFz{37} YWFz{38}
               YWFz{45} YWFz{46} YWFz{47} YWFz{48}]; 

F_matrix1 = inv(Mob_Matrix1);