%% ondas esfericas 3D

clear all;
close all;
clc;

%% inicializacion

%generacion de grilla es importante puesto que en la grilla se pueden
%definir funciones en este caso de dos variables (x,y) es decir que 
%la grilla representa la discretizacion del plano

%empezaremos como base un espacio plano bidimensional, rectangulat de
%dimensiones de Lx = 5(m) Ly = 3(m) con una discretizacion de 0.1(m)

Lx = 20;
Ly = 20;
incx = 0.1;
incy = 0.1;

%generamos la grilla usando el comando meshgrid
[x,y] = meshgrid(0:incx:Lx, 0:incy:Ly);

%las expresiones    0:incx:Lx
%                   0:incx:Lx
%son similares a las usadas en los ciclos for genéricos

%posición de la fuente la fuente está fuera del espacio de cálculo
x0 =  -0.25;
y0 =  10.0;

%distancia a la fuente
r = sqrt( (x-x0).^2 + (y-y0).^2  );

%% incializacion

f = 0.25;               %frecuencia
w = 2*pi*f;             %frecuencia angular
c = 1;                  %velocidad del sonido
k = w/c;                %número de onda
lamda = 2*pi/k;         %longitud de onda
A = 1.5;                %amplitud


%% calculo principal

%inicializamos el tiempo

for t = 0: 0.1 : 2*pi %inicio de ciclo del tiempo (t)
    

    %variable presion sonora p = (r,t)
    p = A*sin(w*t - k*r)./r;  %onda esferica
        
    %graficamos a cada iteración del tiempo se cambia
    %la figura generando la animación
    
    figure(1)
    surf(x,y,p,'FaceColor','interp',...
                'EdgeColor','none',...
                'FaceLighting','phong');
    title('onda esferica')
    xlabel('x (m)');
    ylabel('y (m)');
    zlabel('p(r(x,y),t) (m)');
    grid on;
    box on;
    axis([0 Lx 0 Ly -2*A 2*A]);
    daspect([1 1 1]);
    colormap jet;
    colorbar;
    
    figure(2)
    surf(x,y,p,'FaceColor','interp',...
                'EdgeColor','none',...
                'FaceLighting','phong');
    title('onda esferica')
    xlabel('x (m)');
    ylabel('y (m)');
    zlabel('p(r(x,y),t) (m)');
    grid on;
    box on;
    axis([0 Lx 0 Ly ]);
    daspect([1 1 1]);
    colormap jet;
    colorbar;
    
    pause(0.01);
    
end;                  %fin del ciclo del tiempo