%% prog integral de reyleigh piston rectangular - cuadrado

clearvars; 
close all;
clc;

%este es un programa que determina el campo sonorro de un piston
%rectangular - cuadrado que vibra con velocidad constante normal
%a la superficie

%% datos inciales
f = 5000;               %frecuencia
c = 344;                %velocidad del sonido
w = 2*pi*f;             %frecuencia angular
k = w/c;                %numero de onda
P0 = 2e-5;              %presion de referencia
ro = 1.18;              %densidad del aire

%% tiempo
t = 0;                  %tiempo referencial

%% prámetros de la fuente rectangular - cuadrada
Ly = 5e-1;              %ancho
Lz = 5e-1;              %alto
Ny = 101;               %numero de fuentes en y
Nz = 101;               %numero de fuentes en z
dy = Ly/Ny;             %elemento de ancho
dz = Lz/Nz;             %elemento de alto
dS = dy*dz;             %elemento de superficie
U0 = 1e-3;              %velocidad
dQ = dS*U0;             %elemento de poder de fuente
R = 0.5;                %radio del circulo

%% elemento de presion - amplitud de fuente simple
A = (1i*ro*c*dQ*k)/(2*pi);

%% posicion fuente circular

x0(1:Ny*Nz,1) = -0.01;
y0(1:Ny*Nz,1) = 0;
z0(1:Ny*Nz,1) = 0;

%contador
n = 1; 

for ny = 1:Ny
    for nz = 1:Nz
        y0(n,1) = (ny -1)*dy - Ly/2 + dy/2;
        z0(n,1) = (nz -1)*dy - Lz/2 + dz/2;
        n = n + 1;
    end;
end;

N = length(y0);
x0(1:N,1) = -0.01;

%figura fuente en el plano (y,z)
figure(1)
plot(y0,z0,'b o');
title('geometria de la fuente 2D');
xlabel('y_0(m)');
ylabel('z_0(m)');
legend('(y_0,z_0)');
axis([-1.1*R,1.1*R,-1.1*R,1.1*R]);
daspect([1,1,1]);
grid on;
box on;

%figura fuente en el espacio (x,y,z)
figure(2)
plot3(x0,y0,z0,'b o');
title('geometria de la fuente 3D');
xlabel('x_0(m)');
ylabel('y_0(m)');
zlabel('z_0(m)');
legend('(x_0,y_0,z_0)');
axis([-0.1,0.1,-1.1*R,1.1*R,-1.1*R,1.1*R]);
daspect([1,1,1]);
grid on;
box on;

%% plano de recepcion 1(x,y,0)

%dimensiones e incremento
xmin1 =  0.0;
incx1 =  0.1;
xmax1 = 20.0;

ymin1 = -10.0;
incy1 =   0.1;
ymax1 =  10.0;

%grilla plano (x1,y1,0) y altura z1
[x1,y1] = meshgrid(xmin1:incx1:xmax1,ymin1:incy1:ymax1);  
z1 = 0;

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

%% distancia desde (x0,y0,z0) a (x1,y1,z1=0)

%prelocalizacion
rn1(1:N,1:Nx1,1:Ny1) = 0;

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

%% presion sonora compleja - modulo de la presion sonora - nivel de presion

P1c(1:Nx1,1:Ny1) = 0;       %presion sonora compleja
P1(1:Nx1,1:Ny1) = 0;        %modulo de la presion sonora

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

P1rms = abs(P1c)/sqrt(2);
Lp1 = 20*log10((P1rms+P0)/P0);

%% graficos plano (x1,y1,z1=0)

NC = 100;       %numero de contornos

%presion rms
figure(3)
surf(x1,y1,P1rms,'FaceColor','interp',...
                'EdgeColor','none',...
                'FaceLighting','phong');
title('presion rms plano (x_1,y_1,z_1)');
xlabel('x_1 (m)');
ylabel('y_1 (m)');
legend('Prms_1(x_1,y_1,z_1)');
axis([xmin1,xmax1,ymin1,ymax1]);
daspect([1,1,1]);
colormap jet
colorbar;
grid on;
box on;

%nivel de presion sonora
figure(4)
surf(x1,y1,Lp1,'FaceColor','interp',...
                'EdgeColor','none',...
                'FaceLighting','phong');
title('nivel de presion sonora (x_1,y_1,z_1)');
xlabel('x_1 (m)');
ylabel('y_1 (m)');
legend('L_{p1}(x_1,y_1,z_1)');
axis([xmin1,xmax1,ymin1,ymax1]);
daspect([1,1,1]);
colormap jet
colorbar;
grid on;
box on;

%presion rms axial
figure(5)
plot(x1(round(Nx1/2),:),P1rms(round(Nx1/2),:));
title('presion rms axial (x_1,0,0)');
xlabel('x_1 (m)');
ylabel('Prms(x_1,0,0) (Pa)');
legend('Prms_{ax}');
axis([xmin1,xmax1,min(P1rms(round(Nx1/2),:)),1.1*max(P1rms(round(Nx1/2),:))]);
grid on;
box on;

%nivel de presion axial
figure(6)
plot(x1(round(Nx1/2),:),Lp1(round(Nx1/2),:));
title('nivel de presion axial (x_1,0,0)');
xlabel('x_1 (m)');
ylabel('Lp(x_1,0,0) (dB)');
legend('Prms_{ax}');
axis([xmin1,xmax1,min(Lp1(round(Nx1/2),:)),1.1*max(Lp1(round(Nx1/2),:))]);
grid on;
box on;