%% CGEFEM00120240131 - 2D Acoustics Space - Nuemann Condition 01 BLK

%01 BLK

clear all; 
close all; 
clc;

%% frequency
fini = 1;
fend = 1000;
df = 0.1;
f = transpose(fini:df:fend);
w = 2*pi*f;
NF = length(f);

%% dimension 01

%dimensions
Lx01 = 1.00;
Ly01 = 2.00; 

%% elements 01

%number of elements 01
number_elements_x01 = 200;
number_elements_y01 = 400;
number_total_elements01 = number_elements_x01*number_elements_y01;

%dim elements 01
dx01 = Lx01/number_elements_x01;
dy01 = Ly01/number_elements_y01;
ax01 = dx01/2;
ay01 = dy01/2;

%% nodes 01

%number of nodes 01
number_nodes_elements01 = 4;
number_nodes_x01 = number_elements_x01 + 1;
number_nodes_y01 = number_elements_y01 + 1;
number_total_nodes01 = number_nodes_x01*number_nodes_y01;

number_nodes_elements = number_nodes_elements01;

%dof nodes 01
dof_nodes = 1;

%number of dirichelet conditions 01
number_diri_condt = 0;

%% coordinates nodes for fem

%coordinates nodes for fem
coord_node = 2;

% part 01
xycoord01 = zeros(number_total_nodes01,coord_node);

%generation coord fem part 01
nnod = 1;
for ny = 1:number_nodes_y01
    for nx = 1:number_nodes_x01
        xycoord01(nnod,1) = (nx-1)*dx01;
        xycoord01(nnod,2) = (ny-1)*dy01;
        nnod = nnod + 1;
    end
end

%total coordinates
xycoord = xycoord01;

%% part 1

%total elements
number_total_elements = number_total_elements01;
                       
%total nodes 
number_total_nodes = number_total_nodes01; 

%total dof 
number_total_dof = dof_nodes*number_total_nodes - ...
                   number_diri_condt;

%% id matrix dririchelet /neumann condition

%inicio
id = zeros(number_total_dof,dof_nodes);

%nodes with dirichelet conditions (ninguno en este caso)
nodDr = ones(number_total_dof,dof_nodes);

cont = 1;   %counter

for n = 1:number_total_dof
    if nodDr(n) == 0
        id(n,1) = 0;
    else
        id(n,1) = cont;
        cont = cont + 1;
    end
end    


%% conectivity matrix

%conectivity matrix 01

conec01 = zeros(number_total_elements01,number_nodes_elements01);

%first row conectivity matrix 01
for nel = 1:number_elements_x01
    conec01(nel,1) = nel;                         
    conec01(nel,2) = nel + 1;
    conec01(nel,3) = nel + number_elements_x01 + 2; 
    conec01(nel,4) = nel + number_elements_x01 + 1;
end 

%other rows conectivity matrix 01
for ny = 1:number_elements_y01 - 1
    for nx = 1:number_elements_x01
        nel = ny*number_elements_x01 + nx;
        conec01(nel,1) = conec01(nel-number_elements_x01,4);
        conec01(nel,2) = conec01(nel-number_elements_x01,3);
        conec01(nel,3) = conec01(nel,2) + number_elements_x01 + 1;
        conec01(nel,4) = conec01(nel,1) + number_elements_x01 + 1;
    end
end

conec = conec01;

%% elemental matrices

%data
ax = ax01;
ay = ay01;
number_nodes_elements = number_nodes_elements01;

%wave velocity
c = 340;

%elemental mass matrix
Me = [4, 2, 1, 2;
      2, 4, 2, 1;
      1, 2, 4, 2;
      2, 1, 2, 4];

Me = (ax*ay)/(9*c^2)*Me;

%elemental stifness matrix

Ke = [(  1/(3*ax^2) + 1/(3*ay^2) ),( -1/(3*ax^2) + 1/(6*ay^2) ),( -1/(6*ax^2) - 1/(6*ay^2) ),(  1/(6*ax^2) - 1/(3*ay^2) );
      ( -1/(3*ax^2) + 1/(6*ay^2) ),(  1/(3*ax^2) + 1/(3*ay^2) ),(  1/(6*ax^2) - 1/(3*ay^2) ),( -1/(6*ax^2) - 1/(6*ay^2) );
      ( -1/(6*ax^2) - 1/(6*ay^2) ),(  1/(6*ax^2) - 1/(3*ay^2) ),(  1/(3*ax^2) + 1/(3*ay^2) ),( -1/(3*ax^2) + 1/(6*ay^2) );
      (  1/(6*ax^2) - 1/(3*ay^2) ),( -1/(6*ax^2) - 1/(6*ay^2) ),( -1/(3*ax^2) + 1/(6*ay^2) ),(  1/(3*ax^2) + 1/(3*ay^2) )];

Ke = (ax*ay)*Ke;


%% global matrices

%mass matrix
M = sparse(number_total_dof,number_total_dof);

%stiffness matrix
K = sparse(number_total_dof,number_total_dof);

% %main loop
for nel = 1:number_total_elements
    
    for noi = 1:number_nodes_elements
        for ngi = 1:dof_nodes
             for noj = 1:number_nodes_elements
                 for ngj = 1:dof_nodes
                    
                     NN = id( conec(nel,noi) , ngi );
                     MM = id( conec(nel,noj) , ngj );
                     
                     if (NN ~= 0) && (MM ~= 0) 
                         
                         fila = dof_nodes*(number_nodes_elements-(number_nodes_elements-1))*(noi - 1) + ngi;
                         colu = dof_nodes*(number_nodes_elements-(number_nodes_elements-1))*(noj - 1) + ngj;

                         M(NN,MM) = M(NN,MM) + Me( fila, colu );
                         K(NN,MM) = K(NN,MM) + Ke( fila, colu );
                     
                     end
                     
                 end
             end
        end 
    end    
end    

%proportional damping matrix
a = 1.00e-6;
b = 1.00e-6;

C = a*M + b*K;

%% forced dirichelet condition

%none

%% eigenvalues and eigenvectors

Nhat = 100;
[PHI,LAM] = eigs(K,M,Nhat,150);

%orthonormalization
norm = transpose(PHI)*M*PHI;

for n = 1:Nhat
    PHI(:,n) = PHI(:,n)/sqrt(norm(n,n));
end

XI = transpose(PHI)*C*PHI;

%natural frequencies modal damping

%modal freq
f_FEM = diag( real(sqrt(LAM)))/(2*pi);
w_FEM = 2*pi*f_FEM;

%modal damping
x_FEM = diag(XI)./w_FEM/2;
x_FEM(1) = 0;

%% frequency response

%receptancy_rs
r = 1;
s = 1;
alpha_rs = zeros(NF,1);

for nf = 1:NF
    
    %FRF receptancy
    for n = 2:Nhat
        alpha_rs(nf,1) = alpha_rs(nf,1) + ( PHI(r,n)*PHI(s,n) )/(-w(nf)^2 + 1i*w(nf)*XI(n,n) + LAM(n,n) );
    end    
  
end

%normalizacion
alpha_rs = alpha_rs/max(abs(alpha_rs));

%% figure frf

figure(1)
semilogy( f, abs(alpha_rs),'b','LineWidth',1);
title(['receptancy \alpha(',num2str(r),',',num2str(s),')']);
xlabel('frequency (Hz)');
ylabel('\alpha (Pa s/m)');
legend(['\alpha(',num2str(r),',',num2str(s),')']);
axis([ f(1), f(NF), 0.1*min(abs(alpha_rs)), 10*max(abs(alpha_rs)) ] );
grid on;
box on;

%% impulse response

tini = 0;
tend = 5;
fs = 44100;
dt = 1/fs;
t = transpose(tini:dt:tend);
NT = length(t);

h_rs = zeros(NT,1);

for nt = 1:NT
    
    %IRF receptancy
    for n = 2:Nhat
       
        h_rs(nt,1) = h_rs(nt,1) + ...
                   ( PHI(r,n)*PHI(s,n) )/( w_FEM(n)*sqrt(1 - x_FEM(n)^2) )*exp(-w_FEM(n)*x_FEM(n)*t(nt))*sin(w_FEM(n)*sqrt(1 - x_FEM(n)^2)*t(nt));
    end    
  
end

%normalizacion
h_rs = h_rs/max(abs(h_rs));

%% figure irf

figure(2)
plot(t, h_rs,'b','LineWidth',1);
title(['receptancy h(',num2str(r),',',num2str(s),')']);
xlabel('time (Hz)');
ylabel('h (Pa s/m)');
legend(['h(',num2str(r),',',num2str(s),')']);
axis([ t(1), t(NT), -1.1*max(abs(h_rs)), 1.1*max(abs(h_rs)) ] );
grid on;
box on;

%sound
sound(h_rs,fs);