%% Prog021FEMRecinto2DCerradoParedesRigidasGeometria

clear all;
close all;
clc;

%% geometry

% polygon type 
ptype = 2;
nside = 5;

% plane XY

% coordinates of polygon 
xt(1,1) = -2;
xt(2,1) =  2;
xt(3,1) =  5;
xt(4,1) =  0;
xt(5,1) = -5;


yt(1,1) =  0;
yt(2,1) =  0;
yt(3,1) =  5;
yt(4,1) =  10;
yt(5,1) =  5;

% axis Z
zt =       4;

% geometry
gd = [ptype;nside;xt;yt];

% decomposed geometry matrix
g = decsg(gd);

% boundary conditions - neumann;

xb = [1;0;1;1;48;48;48;48;49;48];

b = [xb xb xb xb xb];  

% plot geometry

figure(101);
pdegplot(g);
axis([-8,8,-2,14]);
grid on;
title('Geometry');
xlabel('x (m)');
ylabel('y (m)');
daspect([1 1 1]);

% mesh
[p,e,t] = initmesh(g);
[p,e,t]=refinemesh(g,p,e,t);
[p,e,t]=refinemesh(g,p,e,t);

figure(102);
pdemesh(p,e,t);
axis([-8,8,-2,14]);
grid on;
title('Mesh');
xlabel('x (m)');
ylabel('y (m)');
daspect([1 1 1]);

% pde coeficients
a = 0;
c = 1;
d = 1;

% eigenproblem plane XY to give natural wave numbers
maxeig = 7;
r = [0 maxeig];
[phi,lamda]=pdeeig(b,p,e,t,c,a,d,r);

% axial and tangencials wave numbers plane XY
kxy = sqrt(lamda);

% matrix dimension
sz = size(phi);
N = sz(1,1);
Nred = sz(1,2);

% matrix dimension
Nfreq = 200;

% material
ro = 1.18;
tau = 1.396364800000000e+005;
Cvel = sqrt(tau/ro);
Radio = 0.03;
S = pi*Radio^2;
P0 = 2e-5;

figure(103);
subplot(2,2,1);
pdesurf(p,t,phi(:,2));
axis([-8,8,-2,14]);
grid on;
title('1st Axial/Tangencial Modal Shape');
xlabel('x (m)');
ylabel('y (m)');
daspect([1 1 1]);
colormap('jet');
box;

subplot(2,2,2);
pdesurf(p,t,phi(:,3));
axis([-8,8,-2,14]);
grid on;
title('2nd Axial/Tangencial Modal Shape');
xlabel('x (m)');
ylabel('y (m)');
daspect([1 1 1]);
colormap('jet');
box;

subplot(2,2,3);
pdesurf(p,t,phi(:,4));
axis([-8,8,-2,14]);
grid on;
title('3th Axial/Tangencial Modal Shape');
xlabel('x (m)');
ylabel('y (m)');
daspect([1 1 1]);
colormap('jet');
box;

subplot(2,2,4);
pdesurf(p,t,phi(:,5));
axis([-8,8,-2,14]);
grid on;
title('4th Axial/Tangencial Modal Shape');
xlabel('x (m)');
ylabel('y (m)');
daspect([1 1 1]);
colormap('jet');
box;

% eigenproblem axis Z
for n1 = 1:Nred
    kz(n1,1) = (n1 -1)*pi/zt;
end;    

% frequency response

% coordinates plane XY
kk = 1;
ll = 3;

% coordinates axis Z
zl = 0;
zk = zt;

%computations

sp(1:Nfreq) = 0;

freq = (1:Nfreq);
w = 2*pi*freq;
        
i1 = ones(1,Nred);
i2 = ones(Nred,1);
    
kxy2 = kxy.^2;
kz2 = transpose(kz.^2);
    
Mkxy2 = kron(kxy2,i1);
Mkz2  = kron(i2,kz2);
    
w0 = Cvel*sqrt(Mkxy2 + Mkz2);
    
BBkl = phi(kk,:).*phi(ll,:);
CCkl = cos(kz*zk).*cos(kz*zl);
Akl = kron(transpose(BBkl),transpose(CCkl));

for n1 = 1:Nfreq
    
    sp(n1) = 0;
    
    for n2 = 1:Nred
        
        for n3 = 1:Nred
            
            sp(n1) = sp(n1) + (j*ro*w(n1)*S*Cvel^2)*(Akl(n2,n3)/(w0(n2,n3)^2 - w(n1)^2));
            
        end;
        
    end;

end;

ff00 = w0/2/pi;

Lp = 20*log10(abs(sp)/P0);

figure(104);
plot(freq,Lp);
axis([20,200,0,140]);
grid on;
title('Sound Pressure Level');
xlabel('f (Hz)');
ylabel('Lp (dB)');
legend('Lp')