%% programa 003 Onda Esferica Modelo Complejo 

clearvars;
close all;
clc;

%graficaremos una onda plana en una caso 2D usaremos una grilla a fin de
%definir las variables (x,y) además del tiempo. 
%la grilla representa una discretizacion del plano

%% parametros principales

%parametros acusticos
c = 344.0;          %velocidad del sonido
f = 1000;           %frecuencia
w = 2*pi*f;         %frecuencia angular
T = 1/f;            %periodo
la = c/f;           %longitud de onda
k = 2*pi/la;        %numero de onda
A = 1.5;            %amplitud
ro = 1.18;          %densidad del aire

%parametros del espacio
Lx = 5;     %limite espacio x
Ly = 5;     %limite espacio y
dx = 0.1;   %incremento en x
dy = 0.1;   %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.5;
y1 =  -0.5;
 
%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.5;
y2 =  0.5;
 
%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;

%posición de la fuente 3 fuera del lugar de estudio
x3 =  5.5;
y3 =  1.0;
 
%programar la distancia a la fuente
%el operador .^ eleva cada elemento en particular de las matrices x e y
r3 = sqrt( (x - x3).^2 + (y - y3).^2 );

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

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(1)
surf(x,y,P,'FaceColor','interp',...
                'EdgeColor','none',...
                'FaceLighting','phong');
title('presion sonora')
xlabel('x(m)');
ylabel('y(m)');
zlabel('P(Pa)');
daspect([1 1 1]);
colormap jet;
grid on;
box on;
colorbar;

%nivel de presion sonora
figure(2)
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 30 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(3)
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(4)
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;