% Program to compute deformation in hypothetical flows 

clear;

% Define the domain, wavenumbers and parameters for defining various
% streamfunction fields

x_domain=3200;
y_domain=4000;

kx=2*pi/x_domain;
ky=2*pi/y_domain;

alpha=ky/2;
beta=0.01;

xc=1000;
yc=y_domain/2;

delta=20;   %grid interval spacing
dplt=10;    %skip interval for plotting vectors

%create the numerical grid

[xp,yp]=meshgrid(0:delta:x_domain-delta,0:delta:y_domain-delta);

% retrieve the grid dimensions
[ny,nx]=size(xp);

% Define perturbation streamfunction

% sinusoidal wave with uniform windspeeds
%psi=-(10^3)*( 6*cos(kx.*(xp-xc)));

% localized straight jet
%psi=-(10^3)*6*sin(2*alpha.*(yp-yc));

% square wave with deformation
%psi=-(10^3)*( 6*cos(kx.*(xp-xc))).*sin(ky.*(yp-yc));

% "realistic case"
psi=-(10^3)*( 6*cos(kx.*(xp-xc-0.2*yp)).*(cos(alpha.*(yp-yc))).^4);

% Define a mean streamfunction corresponding to uniform westerly flow
psi_mean=10^3*(50-beta*yp);
     

%compute velocity fields using "circshift" to deal with the periodic BC
% Note how the indices are shifted to compute x and y derivatives

vs=(circshift(psi,[0,-1])-circshift(psi,[0,+1]))./(2*delta);
us=-(circshift(psi,[-1,0])-circshift(psi,[1,0]))./(2*delta);

% mean shear is uniform, so we just compute a scalar value with a one-sided
% difference
us_mean=-(psi_mean(2,1)-psi_mean(1,1))/delta;

% compute the total wind field
us=us+us_mean;

% Define values of velocity on a more widely spaced grid for plotting
% vectors
[xplt,yplt]=meshgrid(0:delta*dplt:x_domain-delta,0:delta*dplt:y_domain-delta);
usplt=us(1:dplt:ny,1:dplt:nx);
vsplt=vs(1:dplt:ny,1:dplt:nx);

% Plot the stream function and the velocity vectors
figure (1);
[C,h1] = contour(xp, yp, psi+psi_mean);
set(h1,'ShowText','on','TextStep',get(h1,'LevelStep')*2)
hold on;
quiver(xplt,yplt,usplt,vsplt);
axis equal tight;
hold off;

% the following may be useful for your last plot

% plot the stream function and plot the axis of dilatation as a vector
% figure(5)
% contour(xp,yp,psi+psi_mean);
% hold on;
% quiver(xplt, yplt, xdefplt,ydefplt, 'ShowArrowHead','off');
% axis equal tight;
% hold off;




