clear all; % Remove items from workspace
clc; % Clear Command Window
close all; % Close all figures

% u = 2yt
% v = x
% ax = 2y + v 2t = 2y+2xt = 2(xt+y)
% ay = u = 2yt
% vorticity =1/2( dv/dx - du/dy) =1/2 (1 - 2t)


[x,y] = meshgrid(linspace(-5,5,101));
t = 0:0.2:2;

npx=size(x,2);
npy=size(y,1);
nt=length(t);

for k=1:nt
    vx(:,:,k) = 2*y*t(k);
    vy(:,:,k) = x;
    ax(:,:,k) = 2*(y+x*t(k));
    ay(:,:,k) = 2*y*t(k);
    vorticity(:,:,k) = 0.5*(1 - 2*t(k));
end
       
% il ciclo precedente è equivalente a
% for i=1:npx
%     for j=1:npy
%         for k=1:nt
%             vx(i,j,k) = 2*y(j)*t(k);
%             vy(i,j,k) = x(i);
%         end
%     end
% end

figure
for k=1:nt
    C=contourf(x,y,sqrt(vx(:,:,k).^2+vy(:,:,k).^2));
    clabel(C);
    title(sprintf('%0.1f sec',t(k)));
    hold on;
    xy = stream2(x,y,vx(:,:,k),vy(:,:,k),[(-5:1:5)*0-5,(-5:1:5)*0+5,(-5:1:5),(-5:1:5)], [(-5:1:5),(-5:1:5),(-5:1:5)*0-5,(-5:1:5)*0+5]);
    streamline(xy);
    % streamline (per 4,2)
    xyA = stream2(x,y,vx(:,:,k),vy(:,:,k),4,2);
    hline=streamline(xyA);
    set(hline, 'Color', 'r')
    xyB = stream2(x,y,-vx(:,:,k),-vy(:,:,k),4,2);
    hline=streamline(xyB);
    set(hline, 'Color', 'r');
    plot(4,2,'ro')
    pause(1);
    hold off;
end

figure;
s=1:5:101;
for k=1:nt
    quiver(x(s,s),y(s,s),vx(s,s,k),vy(s,s,k),5);
    xlim([-5 5]);
    ylim([-5 5]);
    title('velocity')
    pause(1);
end