%% -------------------------------
%% EPS 131: Intro to Phys. Ocean.
%% Laure Zanna 
%% Matlab Session : Nov. 18th
%% Very useful online help 
%% http://www.mathworks.com/access/helpdesk/help/helpdesk.shtml
%% -------------------------------

clear all; close all; 

%--------------------------------------------------------------------------
% 1. velocity in 2D using quiver(X,Y,u,v)
%--------------------------------------------------------------------------
% grid
x=-1:0.1:1; y=-1:0.1:1;
[X,Y]=meshgrid(x,y);

u=-Y.*exp(-sqrt(X.^2+Y.^2));
v=X.*exp(-sqrt(X.^2+Y.^2));
figure
axes('position', [0.3,0.5,0.35,0.35]); %[left, bottom, width, height])
quiver(X,Y,u,v)
xlabel('x','fontsize',13); 
ylabel('y','fontsize',13); 

%--------------------------------------------------------------------------
% 2. velocity in 3D using quiver3(X,Y,Z,U,V.W)
%--------------------------------------------------------------------------
[X,Y] = meshgrid(-2:0.25:2,-1:0.2:1);
Z = X.* exp(-X.^2 - Y.^2);
[U,V,W] = surfnorm(X,Y,Z);
axes('position', [0.3,0.1,0.35,0.35]); %[left, bottom, width, height])
quiver3(X,Y,Z,U,V,W,0.5);
view(-35,45)
axis ([-2 2 -1 1 -.6 .6])

% ---------------------------------------------------------------
% 3. wind stress with quiver (HW-07)
% ---------------------------------------------------------------
clear all; 
%parameters
tau_0=0.07;
theta=20:4:60; %latitude in deg
length_th=length(theta); %length of the array 
phi=zeros(1,length_th); % define an array of the same length of theta containing only zeros

% grid
[Theta,Phi]=meshgrid(theta,phi);

% wind stress
tau_x=tau_0*cos((pi/20)*(40-Theta));
tau_y=0*Theta;

figure; 
quiver(Phi,Theta,tau_x,tau_y)
set(gca,'XLim',[-4 4],'YLim',[18 62]) %set the limits of the axis
set(gca,'XTick',-4:4:4,'YTick',[20 30 40 50 60]);
set(gca,'XTickLabel',{'-0.8','0','0.8'},'YTicklabel',{'20N','30N','40N','50N','60N'}); 
title('Wind stress \tau^x as fct of latitude','fontsize',13)

%--------------------------------------------------------------------------
% 4. Ekman spiral (HW-06)
%--------------------------------------------------------------------------
clear all; 
% grid in x y z 
xx=zeros(1,11); 
yy=zeros(1,11); % define zeros for the x-y direction
zz=0:-4:-60; %define the depth 
[x,y,z]=meshgrid(xx,yy,zz);

clear xx yy zz
%parameters
f=2*(7.3*10^(-5))*sin(45*pi/180); % coriolis parameter
A_v=10^-2; % friction coefficient
delta=(f/(2*A_v))^(1/2); % 1/Ekman depth
 
% velocity field
u=exp(delta*z).*cos(delta*z); 
v=exp(delta*z).*sin(delta*z);
w=0*z; % define w=0 just for the plot
figure; 
axes('position', [0.3,0.3,0.5,0.5]); %[left, bottom, width, height])
% use quiver3 to draw the Ekman spiral
quiver3(x,y,z,u,v,w); 
xlabel('x','fontsize',13);  
ylabel('y','fontsize',13);  
zlabel('z','fontsize',13); 
view(-35,45) % set  the angle of view (azimuth, vertical)
grid off

%--------------------------------------------------------------------------
% 5. Streamlines for the subtropical gyre (Section-07)
%--------------------------------------------------------------------------
clear all;
%parameters
L_x=5000*10^3 ; 
L_y=3330*10^3 ;  
C=7.3547;
%grid
y=0:L_y/20:L_y; 
x=0:L_x/100:L_x;
[X,Y]=meshgrid(x,y);
% streamlines 
psi=-C*sin(pi*Y/L_y).*(X-L_x);
figure; 
contourf(x,y,psi/10^6,20); 
title('\psi [Sv]'); 
colorbar;
set(gca,'XLim',[0 L_x],'YLim',[0 L_y])
set(gca,'XTick',[0 L_x],'YTick',[0 L_y/2 L_y]);
set(gca,'XTicklabel',{'60W','10W'},'YTicklabel',{'20N','35N','50N'}); 

%--------------------------------------------------------------------------
% 6. contours in 3D : contour3(X,Y,Z)
%--------------------------------------------------------------------------
[X,Y] = meshgrid([-2:.25:2]);
Z = X.*exp(-X.^2-Y.^2);
figure; 
contour3(X,Y,Z,30)
grid off
view(-15,25)

% ------------------------------------------------------------
% 7. plot streamlines from the velocity: streamslice(X,Y,u,v)
% -----------------------------------------------------------
%% streamlines eqn: y=C/x
clear all; 
n=10;m=10;
[X,Y] = meshgrid(-10:10,-10:10);
u=-X; v=Y;
streamslice(X,Y,u,v);
xlabel('X','fontsize',13); ylabel('Y','fontsize',13); 
title('Streamlines','fontsize',14)
