clear all;

% constants
f = 1.25:1.25:2000;
roh = 2.7e3; % Density                           
E = 7.1e10;  % Young's modulus
h = 0.01905; % thickness of the plate                        
nue = 0.33;  % Poisson's ratio
lx = 1.5;                                
ly = 2.12;   
x1 = 1.125;                            
y1 = 1.06;                              
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)); 

Yf1z1zsummeges15 = zz;
Yf1z1zsummeges16 = zz;
Yf1z1zsummeges17 = zz;
Yf1z1zsummeges18 = zz;
Yf1z1zsummeges25 = zz;
Yf1z1zsummeges26 = zz;
Yf1z1zsummeges27 = zz;
Yf1z1zsummeges28 = zz;
Yf1z1zsummeges35 = zz;
Yf1z1zsummeges36 = zz;
Yf1z1zsummeges37 = zz;
Yf1z1zsummeges38 = zz;
Yf1z1zsummeges45 = zz;
Yf1z1zsummeges46 = zz;
Yf1z1zsummeges47 = zz;
Yf1z1zsummeges48 = zz;

Yf1z1zsumme15 = zz;
Yf1z1zsumme16 = zz;
Yf1z1zsumme17 = zz;
Yf1z1zsumme18 = zz;
Yf1z1zsumme25 = zz;
Yf1z1zsumme26 = zz;
Yf1z1zsumme27 = zz;
Yf1z1zsumme28 = zz;
Yf1z1zsumme35 = zz;
Yf1z1zsumme36 = zz;
Yf1z1zsumme37 = zz;
Yf1z1zsumme38 = zz;
Yf1z1zsumme45 = zz;
Yf1z1zsumme46 = zz;
Yf1z1zsumme47 = zz;
Yf1z1zsumme48 = 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; 
    fr = wmnr/(2*pi);
    eta = 0.34876*exp(-fr/161.8177)+0.01809;  % loss factor 
    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./(.5*(roh.*h.*lx.*ly).*(wmn.^2-w.^2));  
       
       Yf1z1zsumme15(z,:) = psi1.*psi5.*rr;
       Yf1z1zsumme16(z,:) = psi1.*psi6.*rr;
       Yf1z1zsumme17(z,:) = psi1.*psi7.*rr;
       Yf1z1zsumme18(z,:) = psi1.*psi8.*rr;
       Yf1z1zsumme25(z,:) = psi2.*psi5.*rr;
       Yf1z1zsumme26(z,:) = psi2.*psi6.*rr;
       Yf1z1zsumme27(z,:) = psi2.*psi7.*rr;
       Yf1z1zsumme28(z,:) = psi2.*psi8.*rr;
       Yf1z1zsumme35(z,:) = psi3.*psi5.*rr;
       Yf1z1zsumme36(z,:) = psi3.*psi6.*rr;
       Yf1z1zsumme37(z,:) = psi3.*psi7.*rr;
       Yf1z1zsumme38(z,:) = psi3.*psi8.*rr;
       Yf1z1zsumme45(z,:) = psi4.*psi5.*rr;
       Yf1z1zsumme46(z,:) = psi4.*psi6.*rr;
       Yf1z1zsumme47(z,:) = psi4.*psi7.*rr;
       Yf1z1zsumme48(z,:) = psi4.*psi8.*rr;
      
       Yf1z1zsummeges15(1,:) = Yf1z1zsumme15(z,:)+Yf1z1zsummeges15(1,:);
       Yf1z1zsummeges16(1,:) = Yf1z1zsumme16(z,:)+Yf1z1zsummeges16(1,:);
       Yf1z1zsummeges17(1,:) = Yf1z1zsumme17(z,:)+Yf1z1zsummeges17(1,:);
       Yf1z1zsummeges18(1,:) = Yf1z1zsumme18(z,:)+Yf1z1zsummeges18(1,:);
       Yf1z1zsummeges25(1,:) = Yf1z1zsumme25(z,:)+Yf1z1zsummeges25(1,:);
       Yf1z1zsummeges26(1,:) = Yf1z1zsumme26(z,:)+Yf1z1zsummeges26(1,:);
       Yf1z1zsummeges27(1,:) = Yf1z1zsumme27(z,:)+Yf1z1zsummeges27(1,:);
       Yf1z1zsummeges28(1,:) = Yf1z1zsumme28(z,:)+Yf1z1zsummeges28(1,:);
       Yf1z1zsummeges35(1,:) = Yf1z1zsumme35(z,:)+Yf1z1zsummeges35(1,:);
       Yf1z1zsummeges36(1,:) = Yf1z1zsumme36(z,:)+Yf1z1zsummeges36(1,:);
       Yf1z1zsummeges37(1,:) = Yf1z1zsumme37(z,:)+Yf1z1zsummeges37(1,:);
       Yf1z1zsummeges38(1,:) = Yf1z1zsumme38(z,:)+Yf1z1zsummeges38(1,:);
       Yf1z1zsummeges45(1,:) = Yf1z1zsumme45(z,:)+Yf1z1zsummeges45(1,:);
       Yf1z1zsummeges46(1,:) = Yf1z1zsumme46(z,:)+Yf1z1zsummeges46(1,:);
       Yf1z1zsummeges47(1,:) = Yf1z1zsumme47(z,:)+Yf1z1zsummeges47(1,:);
       Yf1z1zsummeges48(1,:) = Yf1z1zsumme48(z,:)+Yf1z1zsummeges48(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;

YWFz15 = abs(Yf1z1zsummeges15(1,:).*xx);  % T = theta!
YWFz16 = abs(Yf1z1zsummeges16(1,:).*xx);
YWFz17 = abs(Yf1z1zsummeges17(1,:).*xx);
YWFz18 = abs(Yf1z1zsummeges18(1,:).*xx);
YWFz25 = abs(Yf1z1zsummeges25(1,:).*xx);
YWFz26 = abs(Yf1z1zsummeges26(1,:).*xx);
YWFz27 = abs(Yf1z1zsummeges27(1,:).*xx);
YWFz28 = abs(Yf1z1zsummeges28(1,:).*xx);
YWFz35 = abs(Yf1z1zsummeges35(1,:).*xx);
YWFz36 = abs(Yf1z1zsummeges36(1,:).*xx);
YWFz37 = abs(Yf1z1zsummeges37(1,:).*xx);
YWFz38 = abs(Yf1z1zsummeges38(1,:).*xx);
YWFz45 = abs(Yf1z1zsummeges45(1,:).*xx);
YWFz46 = abs(Yf1z1zsummeges46(1,:).*xx);
YWFz47 = abs(Yf1z1zsummeges47(1,:).*xx);
YWFz48 = abs(Yf1z1zsummeges48(1,:).*xx);

% loglog (f,YWFz15)
% grid on


Mobm = [YWFz15;YWFz16;YWFz17;YWFz18];
Mobm(:,:,2) = [YWFz25; YWFz26; YWFz27; YWFz28];
Mobm(:,:,3) = [YWFz35; YWFz36; YWFz37; YWFz38];               
Mobm(:,:,4) = [YWFz45; YWFz46; YWFz47; YWFz48];

Invmob = zeros(4,1600,4);
v = ones(4,1600,1);
F = zeros(4,1600,1);
for i = 1:1600 
    Invmob(:,i,:) = pinv(squeeze(Mobm(:,i,:)));
    F(:,i,:) = squeeze(Invmob(:,i,:))*squeeze(v(:,i,:));
end

           
% Y15=YWFz15(:,1);     
% Y16=YWFz16(:,1);
% Y17=YWFz17(:,1);
% Y18=YWFz18(:,1);
% Y25=YWFz25(:,1);
% Y26=YWFz26(:,1);
% Y27=YWFz27(:,1);
% Y28=YWFz28(:,1);
% Y35=YWFz35(:,1);
% Y36=YWFz36(:,1);
% Y37=YWFz37(:,1);
% Y38=YWFz38(:,1);
% Y45=YWFz45(:,1);
% Y46=YWFz46(:,1);
% Y47=YWFz47(:,1);
% Y48=YWFz48(:,1);
% 
% mobm = [Y15 Y16 Y17 Y18;
%         Y25 Y26 Y27 Y28;
%         Y35 Y36 Y37 Y38;
%         Y45 Y46 Y47 Y48];
% 
%     v = [1 1 1 1]; 
% F = v/mobm;


% Mobm = [YWFz15;YWFz16;YWFz17;YWFz18YWFz25;YWFz26;YWFz27;YWFz28;YWFz35;YWFz36;YWFz37;YWFz38;YWFz45;YWFz46;YWFz47;YWFz48];
% Mobm(:,:,2) = [YWFz15; YWFz16];
% loglog(Mobm(1,:,2))



% YWFz17 YWFz18
%         YWFz25 YWFz26 YWFz27 YWFz28         
%         YWFz35 YWFz36 YWFz37 YWFz38        
%         YWFz45 YWFz46 YWFz47 YWFz48];
% 
% 
% v1 = ones(1,3200);
% 
% F1 = v1./YWFz15+v1./YWFz16+v1./YWFz17+v1./YWFz18;
% F2 = v1./YWFz25+v1./YWFz26+v1./YWFz27+v1./YWFz28;
% F3 = v1./YWFz35+v1./YWFz36+v1./YWFz37+v1./YWFz38;
% F4 = v1./YWFz45+v1./YWFz46+v1./YWFz47+v1./YWFz48;