% Visualizacion de un campo de velocidades 2D
% no uniforme y no estacionario. 
% 
% Problema explorado en clase practica (21-Ago-2014)
% Estructura de la Materia 1 - 2do Cuat. 2014
%
% Nota:
% Este codigo tiene como objetivo ilustrar un ejemplo de lineas de 
%  corriente, trayectorias y de trazas de un flujo no uniforme y no 
%  estacionario. El mismo ha sido escrito con la intencion de ser facil-
%  mente legible por sus potenciales usuarios, por lo que no se ha buscado
%  optimizar su tiempo de ejecucion explotando las capacidades de calculo
%  'vectorizado' que ofrece MATLAB.
%
% Pablo Cobelli
% cobelli at df dot uba dot ar
%
% Licencia: Creative Commons Attribution-NonCommercial-ShareAlike 3.0.


% borramos toda variable y cerramos toda figura preexistentes
clear all
close all

% definimos los vectores de espacio 2D y tiempo
x = linspace(0,10,30);
y = linspace(-3,3,25);
t = linspace(0,20,101);

% creamos una grilla cartesiana 2D a partir de los vectores x e y (1D)
[xg, yg] = meshgrid(x,y);

% definimos tambien vectores espaciales con mayor resolucion
%  (para obtener lineas de trayectoria y de trazas mejor resueltas 
%  espacialmente)
xbd = linspace(0,10,200);
ybd = linspace(-3,3,150);

% fijamos valores concretos a los parametros del problema
uo = 0.3;
vo = 0.1;
w  = 0.5;

% fijamos las coordenadas de un punto en el espacio;
%  esto lo usaremos como inicio de trayectorias y tambien como
%  punto por donde pasan las lineas de corriente y de trazas.
%   Atencion: esto es por comodidad!, podriamos utilizar distintos puntos
%   para las distintas formas de visualizar el flujo.
xo = 0;
yo = 0;

% para cada instante de tiempo, hacemos:
for ii=1:numel(t);

    % determinamos el campo de velocidades 2D
    u = uo*ones(size(xg));
    v = vo*cos(w*(t(ii)-xg/uo));
    
    % seteamos la figura #1 como activa
    figure(1);
    
    % graficamos en negro el campo de velocidades 
    quiver(x,y,u,v,0.75,'k');
    
    % evitamos la que proxima funcion de graficacion nos borre el
    % grafico que acabamos de hacer, ya que queremos superponerlos todos.
    hold on;
    
    % graficamos en rojo una linea de corriente [que pasa por (xo, yo)]
    lincorr = streamline(xg,yg,u,v,xo,yo);
    set(lincorr,'Color','red','Linewidth',2)
    
    % graficamos en verde la trayectoria de una particula que a tiempo
    % t = 0 estaba en (xo,yo)
    xtray1 = uo*t(ii);
    ytray1 = 0;
    plot(xtray1,ytray1,'go','MarkerFaceColor','g');
    
    % graficamos en magenta la trayectoria de una particula que a tiempo
    % t = pi/(2*w) estaba en (xo,yo)
    xtray2 = uo*(t(ii)-pi/2/w);
    ytray2 = vo*(t(ii)-pi/2/w);
    plot(xtray2,ytray2,'mo','MarkerFaceColor','m');
    plot(xbd,vo/uo*xbd,'m','Linewidth',2);
    
    % graficamos en azul la lineas de trazas que pasa por (xo,yo)
    yzbd = yo + vo*sin(w*t(ii) - w*xbd/uo).*(xbd-xo)/uo ;
    plot(xbd,yzbd,'b-','Linewidth',2);
    
    % hacemos que el grafico tenga relacion de aspecto 1:1
    daspect([1 1 1]);
    
    % le agregamos una grilla para referencia
    grid on;
    
    % fijamos la extension de los ejes en ambos sentidos 
    % para no desperdiciar espacio en la figura
    axis([min(x) max(x) min(y) max(y)]);
    
    % le agregamos una leyenda al grafico
    legend('u(x,y,t)','Linea de corriente por (0,0)','Particula 1',...
        'Particula 2','Trayectoria particula 2','Traza por (0,0)',...
        'location','southwest');    
    
    % agregamos titulo al grafico (con un rotulo temporal),
    % y tambien titulos a los ejes de coordenadas
    title(['T = ' num2str(t(ii)) ' sec']);
    xlabel('x [m]');
    ylabel('y [m]');
    
    % exigimos que MATLAB actualice la imagen en pantalla
    % para asegurarnos que grafico lo que le pedimos hasta aqui
    drawnow;
    
    % liberamos el grafico, de forma tal de que todo lo graficado
    % de ahora en mas no se superpone a lo anterior, sino que el grafico
    % anterior se reemplaza por otro nuevo
    hold off;
    
    % introducimos una pausa de 0.1 sec para tener tiempo de ver el 
    % grafico en pantalla antes de que se pase al siguiente instante
    % de tiempo
    pause(0.1);
    
end