%% visualization of curl and divergence

[x,y] = meshgrid(-10:1:10);  % for vector field
[x2,y2] = meshgrid(-10:0.2:10);  % for images

map = [[linspace(0,1,128)';ones(128,1)],[linspace(0,1,128)';linspace(1,0,128)'],[ones(128,1);linspace(1,0,128)']];

%% vector fields from sheet 1

% curl-free, not divergence-free

vx = x;
vy = y;

subplot(121)
quiver(x,y,vx,vy)
xlabel('x')
ylabel('y')
axis image
title('f = (x,y) : curl f = 0, div f = 2')

% divergence-free, not curl-free

vx = -y;
vy = x;

subplot(122)
quiver(x,y,vx,vy)
xlabel('x')
ylabel('y')
axis image
title('f = (-y,x) : curl f = 2, div f = 0')

%% function with both curl and divergence

vx = x.^2 - y.^2;
vy = 2*x.*y;

div = 4*x2;
curl = 4*y2;

figure

subplot(121)
hold on
imagesc(-10:10,-10:10,div,[-40 40])
quiver(x,y,vx,vy,'k')
xlabel('x')
ylabel('y')
axis image
hold off
title('div f = 4x')
colorbar

subplot(122)
hold on
imagesc(-10:10,-10:10,curl,[-40 40])
quiver(x,y,vx,vy,'k')
xlabel('x')
ylabel('y')
axis image
hold off
title('curl f = 4y')
colorbar

% hack for global title (since suptitle is only available in some toolbox)
a = axes;
t1 = title('f(x,y) = (x^2-y^2, 2xy)');
a.Visible = 'off';
t1.Visible = 'on';

colormap(map)

%% more complicated function with both curl and divergence

vx = cos(0.2*x+0.4*y);
vy = sin(0.2*x-0.4*y);

div = -.2*sin(.2*x2+.4*y2)-.4*cos(.2*x2-.4*y2);
curl = .2*cos(.2*x2-.4*y2)+.4*sin(.2*x2+.4*y2);

figure

subplot(121)
hold on
imagesc(-10:10,-10:10,div)
quiver(x,y,vx,vy,'k')
xlabel('x')
ylabel('y')
axis image
% scatter(5,-5,1000,'w')
hold off
title('div f')
colorbar

subplot(122)
hold on
imagesc(-10:10,-10:10,curl)
quiver(x,y,vx,vy,'k')
xlabel('x')
ylabel('y')
axis image
% scatter(4,1.25,1000,'w')
hold off
title('curl f')
colorbar

% hack for global title (since suptitle is only available in some toolbox)
a = axes;
t1 = title('f(x,y) = (cos(0.2x+0.4y), sin(0.2x-0.4y))');
a.Visible = 'off';
t1.Visible = 'on';

colormap(map)

%%

clear x y vx vy div curl