%% prog integral de Reyleygh piston cuadrado

clearvars;
close all;
clc;

%este programa calcula a partir de una aproximación de la integral de
%reyleigh el campo sonoro producido por un pistón cuadrado/rectangular
%que vibra con una velocidad constante

%% datos inciales

f = 5000;       %frecuencia
c = 344;        %velocidad del sonido
w = 2*pi*f;     %frecuencia angular
la = c/f;       %longitud de onda
k = w/c;        %numero de onda
P0 = 2e-5;      %presion de referencia
ro = 1.18;      %densidad del aire
%% tiempo
t = 0;          %tiempo referencial

%% parametros fuente cuadrada/rectangular

%dimensiones
Ly = 5e-1;      %ancho
Lz = 5e-1;      %alto
Nfy = 51;       %numero de fuentes en y
Nfz = 51;       %numero de fuentes en z
dy = Ly/Nfy;    %elemento de longitud en y
dz = Lz/Nfy;    %elemento de longitud en z
dS = dy*dz;     %elemento de superficie
U0 = 5e-3;      %velocidad constante en todo el piston
dQ = dS*U0;     %elemento de poder de fuente

%amplitud de fuente simple
A = (j*ro*c*dQ*k)/(2*pi);

%% posicionamiento de fuentes simples que forman el piston cuadrado/rectangular

%prelocalizacion
x0(1:Nfy*Nfz) = -0.01;
y0(1:Nfy*Nfz) =  0;
z0(1:Nfy*Nfz) =  0;

%contador
n = 1;

%posicionamiento de fuentes simples
for ny = 1:Nfy
    for nz = 1:Nfz
        y0(n) = (ny-1)*dy - Ly/2 + dy/2;
        z0(n) = (nz-1)*dz - Lz/2 + dz/2;
        n = n + 1;
    end;
end;

N = length(y0);

figure(1)
plot(y0,z0,'o');
title('geometria de la fuente');
xlabel('y0 (m)');
xlabel('z0 (m)');
axis([-1.5*Ly/2,1.5*Ly/2,-1.5*Lz/2,1.5*Lz/2]);
daspect([1 1 1]);
grid on;
box on;

%% plano de recepcion (X,Y,0)        

%dimensiones
xmin = 0;
xmax = 20;
incx = 0.05;

ymin = -10;
ymax =  10;
incy = 0.05;

%grilla plano xy
[x1,y1] = meshgrid(xmin:incx:xmax,ymin:incy:ymax);
z1 = 0;                 %altura plano receptor

%longitud de la grilla
Nx1 = length(x1);       %longitud de la grilla en x1
Ny1 = length(y1);       %longitud de la grilla en y1

%prelocalizacion de las distancias como un arreglo 3D
rn1(1:N,1:Nx1,1:Ny1) = 0; 

%distancias
for n1 = 1:N
    rn1(n1,:,:) = sqrt( (x1-x0(n1)).^2 + (y1-y0(n1)).^2 + (z1-z0(n1)).^2 );
end;

%presion sonora prelocalizacion
P1c(1:Nx1,1:Ny1) = 0; %presion sonora compleja
P1(1:Nx1,1:Ny1) = 0;  %presion sonora en módulo del valor complejo amplitud

for n1 = 1:N
        pn(:,:) =  A*exp(j*(w*t - k*rn1(n1,:,:) ) )./(rn1(n1,:,:));
        P1c = P1c(:,:) + pn;
end;
P1 = abs(P1c);

%nivel de presion sonora
Lp1 = 20*log10(P1/P0/sqrt(2));

%% graficos plano (X,Y,0)

figure(2)
surf(x1,y1,P1,'FaceColor','interp',...
                'EdgeColor','none',...
                'FaceLighting','phong');
title('presion sonora plano xy z = 0');
xlabel('x (m)');
ylabel('y (m)');
zlabel('p (N/m^2)')
colormap jet
colorbar
grid on;
box on;

figure(3)
surf(x1,y1,Lp1,'FaceColor','interp',...
                'EdgeColor','none',...
                'FaceLighting','phong');
title('nivel de presion sonora plano xy z = 0');
xlabel('x (m)');
ylabel('y (m)');
zlabel('Lp (dB)')
colormap jet
colorbar
grid on;
box on;
  
figure(4)
plot(x1( round(Nx1/2),:),P1( round(Nx1/2),:) );
title('presión sonora axial');
xlabel('x (m)');
ylabel('Pax (N/m^2)');
grid on;
box on;
  
figure(5)
plot(x1( round(Nx1/2),:),Lp1( round(Nx1/2),:) );
title('nivel de presión sonora axial');
xlabel('x (m)');
ylabel('Pax (N/m^2)');
grid on;
box on;