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

% u = 2y
% v = -xy
% w = zx;
% ax = -2xy
% ay = -2 y^2 + x^2 y
% az = 2 y z + x^2 z
% vorticity_x = 0 
% vorticity_y = -1/2 z
% vorticity_z = 1/2( dv/dx - du/dy) = -1/2 y -1

velx =@(x,y,z) 2*y;
vely =@(x,y,z) -x.*y;
velz =@(x,y,z) x.*z;
vorx =@(x,y,z) zeros(size(x,3),size(y,2),size(z,1));
vory =@(x,y,z) -0.5*z;
vorz =@(x,y,z) -0.5*y-1;


[x,y,z] = meshgrid(linspace(-5,5,101));

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

Vx = velx(x,y,z);
Vy = vely(x,y,z);
Vz = velz(x,y,z);
ax(1:npx,1:npy,1:npz) = -2*x.*y;
ay(1:npx,1:npy,1:npz) = -2*y.^2+x.^2.*y;
az(1:npx,1:npy,1:npz) = -2*y.*z+x.^2.*z;
vorticity_x = vorx(x,y,z);
vorticity_y = vory(x,y,z);
vorticity_z = vorz(x,y,z);

[s_x, s_y] = meshgrid(linspace(-5,5,11),linspace(-5,5,11));
s_z = ones(11,11);

starting_x =repmat(reshape(s_x,1,11*11),1,3);
starting_y =repmat(reshape(s_y,1,11*11),1,3);
starting_z =[-5*reshape(s_z,1,11*11)  0*reshape(s_z,1,11*11) 5*reshape(s_z,1,11*11)];


figure;
hold on
s=1:20:101;
view(15,30)
quiver3(x(s,s,s),y(s,s,s),z(s,s,s),Vx(s,s,s),Vy(s,s,s),Vz(s,s,s),2);
quiver3(x(s,s,s),y(s,s,s),z(s,s,s),vorticity_x(s,s,s),vorticity_y(s,s,s),vorticity_z(s,s,s),2);
xlim([-5 5]);
ylim([-5 5]);
zlim([-5 5])
title('velocity')
grid on


figure;
hold on;
view(15,30)
h = slice(x,y,z,sqrt(Vx.^2+Vy.^2+Vz.^2),[],[],[-5,0,5]);
set(h,'EdgeColor','none','FaceAlpha',0.75');
h = streamline(x,y,z,Vx,Vy,Vz,starting_x,starting_y,starting_z);
set(h,'Color','r');
h=quiver3(starting_x,starting_y,starting_z,...
    velx(starting_x,starting_y,starting_z),...
    vely(starting_x,starting_y,starting_z),...
    velz(starting_x,starting_y,starting_z),2);
set(h,'Color','b');
h = streamline(x,y,z,vorticity_x,vorticity_y,vorticity_z,starting_x,starting_y,starting_z);
set(h,'Color','g');
h=quiver3(starting_x,starting_y,starting_z,...
    vorx(starting_x,starting_y,starting_z),...
    vory(starting_x,starting_y,starting_z),...
    vorz(starting_x,starting_y,starting_z),2);
set(h,'Color','y');
xlim([-5 5]);
ylim([-5 5]);
zlim([-5 5])
title('velocity')
grid on