% compute u,v using finite differencing with hypothetically 
% observed pressure data and compare with actual values
close all; clear all;
% path to colorbarf:
cygwin_root='/cygwin';  % set to '/cygwin' or ''
addpath([cygwin_root '/home/eli/matlab'],'-end');

dx=0.05;dy=0.05; %spacing between gridpoints
% centers of high/ low:
x1=0.5; y1=0.0;
x2=-0.5; y2=0.0;
% scale of high/ lows:
Lx=0.4; Ly=0.4;
[x y]=meshgrid(-1:dx:1,-1:dy:1); %grid of measuring instruments
% pressure field:
p=exp(-(x-x1).^2/Lx^2-(y-y1).^2/Ly^2) ...
  -exp(-(x-x2).^2/Lx^2-(y-y2).^2/Ly^2); 
[nx,ny]=size(p); %number of x,y gridpoints (this stores to both nx and ny)

%u=-dp/dy, v=dp/dx

%create finite difference computed velocity variables
u_fd=zeros(nx,ny); v_fd=u_fd; 

%loop through all y to compute u=-dp/dy
for n=2:nx-2
    u_fd(n,:)=-(p(n+1,:)-p(n-1,:))/(2*dy); %compute -dp/dy for yn and all x
end

%loop through all x to compute v=dp/dx
for n=2:nx-2
    v_fd(:,n)=(p(:,n+1)-p(:,n-1))/(2*dx); %compute dp/dx for xn and all y
end

figure
s=2; % subsampling of quiver arrows.
[c,h]=contourf(x,y,p);
hold on
colorbarf(c,h);
quiver(x(1:s:nx,1:s:nx),y(1:s:nx,1:s:nx),u_fd(1:s:nx,1:s:nx),v_fd(1:s:nx,1:s:nx));
title('pressure contours and geostrophic velocity arrows')
%clabel(c,h);
xlabel('x (east-west distance)')
ylabel('y (north-south distance)')
axis square

% u=-2*y;v=2*x;
% figure
% quiver(x,y,u,v),title('pressure contours and geostrophic velocity arrows')
% hold on
% contour(x,y,p)
% xlabel('x (east-west distance)')
% ylabel('y (north-south distance)')
