%% Numerical Modelling in Geosciences
%  by Manuele Faccenda, Università di Padova
%
%  Practice 2
%
%  Exercise 1.2 from book "Introduction to numerical geodynamic modelling"
%  by Taras Gerya, ETH Zuerich
 
clear all;
close all
clf;

%% Set,W,H,vx0,vy0
W=1e+6;
H=1.5e+6;
vx0=1e-9;
vy0=vx0;

%2D Grid
xnum = 31;
ynum = 31;
xstp=W/(xnum-1);
ystp=H/(ynum-1);
ax=0:xstp:W;
ay=0:ystp:H;
X=zeros(xnum,ynum);
Y=zeros(xnum,ynum);
for i = 1 : xnum
    for j = 1 : ynum
        if i > 1
            X(i,j)=X(i-1,j)+xstp;
        end
        if j > 1
            Y(i,j)=Y(i,j-1)+ystp;
        end
    end
end

%Initialize vx, vy
vx=zeros(xnum,ynum);
vy=zeros(xnum,ynum);
%Compute vx,vy at every grid point
for i = 1 : xnum
    for j = 1 : ynum
        vx(i,j) = -vx0*sin(2*pi*X(i,j)/W)*cos(pi*Y(i,j)/H);
        vy(i,j) =  vy0*cos(2*pi*X(i,j)/W)*sin(pi*Y(i,j)/H);
    end
end

figure(1)
subplot(1,2,1)
pcolor(ax,ay,rot90(vx))
%surf(ax,ay,vx)
shading interp
colorbar
axis  image
title('Vx')
xlabel('W')
ylabel('H')

subplot(1,2,2)
pcolor(X,Y,vy)
%surf(X,Y,vy)
shading interp
colorbar
axis ij image
title('Vy')
xlabel('W')
ylabel('H')

%% Compute dvx/dx, dvy/dy, and div(v) analytically 
dvxdx=zeros(xnum,ynum);
dvydy=zeros(xnum,ynum);
for i = 1 : xnum
    for j = 1 : ynum
        dvxdx(i,j) = -vx0*cos(2*pi*X(i,j)/W)*2*pi/W*cos(pi*Y(i,j)/H);
        dvydy(i,j) = vy0*cos(2*pi*X(i,j)/W)*cos(pi*Y(i,j)/H)*pi/H;
    end
end

figure(2)
subplot(1,2,1)
pcolor(X,Y,dvxdx)
%surf(X,Y,dvx)
shading interp
colorbar
axis ij image
title('dvx')
xlabel('W')
ylabel('H')

subplot(1,2,2)
pcolor(X,Y,dvydy)
%surf(X,Y,dvy)
shading interp
colorbar
axis ij image
title('dvy')
xlabel('W')
ylabel('H')

%Compute divergence
div = dvxdx + dvydy;
figure(3)
pcolor(X,Y,div)
shading interp
colorbar
axis ij image
title('div')
xlabel('W')
ylabel('H')
hold on
quiver(X,Y,vx,vy)
axis ij image


%Find vx0 and vy0 such that div(v) = 0

% div(v) = dvx/dx + dvy/dy = 0

% -vx0*cos(2*pi*x/W)*2*pi/W*cos(pi*y/H) + vy0*cos(2*pi*x/W)*cos(pi*y/H)*pi/H = 0
%
% vy0*cos(2*pi*x/W)*cos(pi*y/H)*pi/H = vx0*cos(2*pi*x/W)*2*pi/W*cos(pi*y/H)
%
% vy0 = vx0*2*H/W
%
% If H = 1500 km and W = 1000 km --> set vy0 = 3*vx0