%% dipolo acústico - factor direccional

clearvars;
close all;
clc;

%datos iniciales
f = 10000;
w = 2*pi*f;
c = 344;
k = w/c;
d = 0.1;
A = 0.5;
ro = 1.18;

%angulos
th = 0:2*pi/360:2*pi;

%factor direccional
H = abs( sin(k*d/2*sin(th)) );

%patron polar
figure(1)
polar(th,H);
title('patron polar dipolo acustico');
legend([ 'f = ',num2str(f),'(Hz)',],'Location','southoutside');
grid on;
box on;


%parametros del espacio
Lx = 5;     %limite espacio x
Ly = 5;     %limite espacio y
dx = 0.05;   %incremento en x
dy = 0.05;   %incremento en y

%grilla desde -Lx a Lx vs -Ly a Ly
[x,y] = meshgrid(-Lx:dx:Lx , -Ly:dy:Ly);

%% graficaremos la presion sonora en una situación estatica del tiempo esferica

%tiempo y presion sonora
t0 = 0;  %tiempo fijo

%posición de la fuente 1 fuera del lugar de estudio
x1 =  5.25;
y1 =  -0.05;
 
%programar la distancia a la fuente
%el operador .^ eleva cada elemento en particular de las matrices x e y
r1 = sqrt( (x - x1).^2 + (y - y1).^2 );

%presion sonora para tiempo fijo usando el modelo de exponencial compleja
%el operador ./ divide a las matrices elemento a elemento
p1 = A*exp( j*(w*t0 - k*r1) )./r1;


%posición de la fuente 2 fuera del lugar de estudio
x2 =  5.25;
y2 =  0.05;
 
%programar la distancia a la fuente
%el operador .^ eleva cada elemento en particular de las matrices x e y
r2 = sqrt( (x - x2).^2 + (y - y2).^2 );

%presion sonora para tiempo fijo usando el modelo de exponencial compleja
%el operador ./ divide a las matrices elemento a elemento
p2 = -A*exp( j*(w*t0 - k*r2) )./r2;


%niveles
P0 = 2e-5;
P = abs(p1 + p2);
Prms = P/sqrt(2);
Lp = 20*log10(Prms/P0);

%% velocidad de partículas

%determinación del vector de velocidad de particulas a partir del 
%gradiente de la presion
%calcula el gradiente de la funcion p(xy) de forma numérica
%             [ dp(xy)   dp(xy) ]
%nabla(pxy) = | ------   ------ |
%             [   dx       dy   ]

[Ux,Uy] = gradient(P,dx,dy);        %calcula el vector gradiente numerico para cada punto
Ux = -Ux/(ro*w);                    %escala la velocidad de partículas Ux
Uy = -Uy/(ro*w);                    %escala la velocidad de partículas Uy
UxUy = -sqrt( Ux.^2 + Uy.^2 );      %longitud del vector gradiente para cada punto     

%% graficos de presion sonora

%presion sonora
figure(2)
surf(x,y,P,'FaceColor','interp',...
                'EdgeColor','none',...
                'FaceLighting','phong');
title('presion sonora')
xlabel('x(m)');
ylabel('y(m)');
zlabel('P(Pa)');
axis([-Lx Lx -Ly Ly ]);
daspect([1 1 1]);
colormap jet;
grid on;
box on;
colorbar;

%nivel de presion sonora
figure(3)
surf(x,y,Lp,'FaceColor','interp',...
                'EdgeColor','none',...
                'FaceLighting','phong');
title('nivel de presion sonora')
xlabel('x(m)');
ylabel('y(m)');
zlabel('Lp(dB)');
axis([-Lx Lx -Ly Ly 0 120]);
colormap jet;
grid on;
box on;
colorbar;

%% grafico presion sonora contorno 2D con flechas velocidad de particulas 2D

NC = 100;                   %numero de contornos

figure(4)
contour(x,y,P,NC);          %contorno 2D
hold on;                    %se usa para no borrar lo dibujado
quiver(x,y,Ux,Uy);          %dibuja flechas 2D 
title('presion sonora contorno  2D con flechas velocidad de particulas 2D');
xlabel('x (m)');
ylabel('y (m)');
zlabel('P(x,y) (Pa)');
daspect([1 1 1]);
grid on;
box on;
colormap jet;
colorbar;

%% grafico presion sonora contorno 3D con flechas velocidad de particulas 3D

figure(5)
contour3(x,y,P,NC);         %contorno 3D
hold on;                    %se usa para no borrar lo dibujado
quiver3(x,y,P,Ux,Uy,UxUy);	%dibuja flechas 3D
title('presion sonora contorno 3D con flechas velocidad de particulas 3D');
xlabel('x (m)');
ylabel('y (m)');
zlabel('P(x,y) (Pa)');
daspect([1 1 1]);
grid on;
box on;
colormap jet;
colorbar;