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

% u = 2yt
% v = x
% w = 2;
% ax = 2y + v 2t = 2y+2xt = 2(xt+y)
% ay = u = 2yt
% az = 0
% vorticity_x = 0 
% vorticity_y = 0
% vorticity_z = 1/2( dv/dx - du/dy) =1/2 (1 - 2t)


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

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

for k=1:nt
    Vx(:,:,:,k) = 2*y*t(k);
    Vy(:,:,:,k) = x;
    Vz(:,:,:,k) = 2*ones(npx,npy,npz);
    ax(:,:,:,k) = 2*(y+x*t(k));
    ay(:,:,:,k) = 2*y*t(k);
    az(:,:,:,k) = zeros(npx,npy,npz);
    vorticity_x(:,:,:,k) = zeros(npx,npy,npz);
    vorticity_y(:,:,:,k) = zeros(npx,npy,npz);
    vorticity_z(:,:,:,k) = 0.5*(1 - 2*t(k));
end

% starting_z = zeros(5,5); % starting points for streamlines
% [starting_y, starting_z] = meshgrid(...
%     linspace(-0.005,.005,5), ...
%     linspace(-0.005,.005,5));
% figure(1)
 
starting_x = repmat([(-5:1:5)*0-5,(-5:1:5)*0+5,(-5:1:5),(-5:1:5)],1,3);
starting_y = repmat([(-5:1:5),(-5:1:5),(-5:1:5)*0-5,(-5:1:5)*0+5],1,3);
starting_z = [-3*ones(1,44), zeros(1,44), 3*ones(1,44)];


figure;
s=1:20:101;
for k=1:nt
    view(45,120)
    quiver3(x(s,s,s),y(s,s,s),z(s,s,s),Vx(s,s,s,k),Vy(s,s,s,k),Vz(s,s,s,k),1);
    xlim([-5 5]);
    ylim([-5 5]);
    zlim([-5 7])
    title(sprintf('velocity %0.1f sec',t(k)))
    grid on
    pause(1);
end

figure;
for k=1:nt
    view(45,120)
    h = slice(x,y,z,sqrt(Vx(:,:,:,k).^2+Vy(:,:,:,k).^2+Vz(:,:,:,k).^2),[],[],[-3,0,3]);
    set(h,'EdgeColor','none','FaceAlpha',0.75');
    h = streamline(x,y,z,Vx(:,:,:,k),Vy(:,:,:,k),Vz(:,:,:,k),starting_x,starting_y,starting_z);
    xlim([-5 5]);
    ylim([-5 5]);
    zlim([-5 5])
    title(sprintf('velocity %0.1f sec',t(k)))
    grid on
    pause(1);
end