%  2D FLOW ROUND A CIRCULAR CYLINDER THAT FILLS THE FULL DEPTH OF THE 
%FLUID ... WITH AND WITHOUT ROTATION, TO SHOW THE EXACT PRESSURE FIELD
%AND ITS VARIATION ALONG STREAMLINES.
 %                                              P.B. Rhines ii 2004
% flow is same in all cases but pressure is not
% The streamfunction psi(x,y) is the sum of a constant mean velocity
% (U,0) and a dipole mass source-sink pair (the 2d term), which mimics 
% the effect of the cylinder (and you can verify that psi = const. on r=a).
% Laplace's equation del squared psi = 0 has these 'singularity' solutions
% like the field of a point charge or plus/minus charge pair (a dipole)
% in electrostatics; here in a fluid the singularities are point sources or
% sinks of flow, or point vortices.  
%    Once you make the perspective plot use the mouse controls to swivel it
%    around the viewing angle. If it's too slow reduce the grid size ..
%    it's pretty big right not.  Straight contour plots using contour( )
%    or using pcolor( );  shading interp  are also good. Simple line plots
%    plot(p); plot(y,p(50,:)  are good; contoured data misses out a lot
%    of important structure. Note the need to use ./ and .* operators with
%    arrays.


[x y]=meshgrid(-50:50,-50:50);
r=sqrt(x.*x+y.*y);
theta=atan2(y,x);   % this works in 4 quadrants whereas atan does not.
a=15; U=2; f=0.1;
Ro=U/(f*a)
p=0.5*U*U*((2*a*a./(r.*r)).*cos(2*theta)-(a./r).^4)-f*U*(r-a*a./r).*sin(theta);
psi=-U*sin(theta).*(r-a*a./r);
jj=find(r<a);
p(jj)=0.; psi(jj)=0.;
hold off
contour(x,y,psi,25,'k')
hold on
contour(x,y,p,25)
axis equal
figure
meshc(p)
title('2D flow past circular cylinder (extending through full depth), p and {\psi} (black), Ro=1.33')


%  STEADY 2D ROTATING FLOW OVER A MOUNTAIN (A SHORT CYLINDER), UNIFORM
%  DENSITY FLUID
%  L=box-size, a=mountain radius,d = mountain height, H=mean fluid depth
%  g=gravity, lambsq= square of Rossby radius, psi = streamfunction
%  ngrid=number of grid points in either direction.
%
% The two solutions inside and outside r=a must be 
% matched: at r=a make eta continuos and d eta/dr continuous
% This is the smoothest possible solution and if we did not do this
% there would be infinite vorticity at r=a; we know from PV conservation 
% that the vorticity is finite  Eqn of motion is
%    delsq psi = fh'/H^2  so zeta definitely has jump across r=a but 
%    azimuthal velocity is continuous. 
%  Solution is eta= A+B r^2  r < a
%                 = C ln(r)  r >= a
%  to which we add a uniform zonal flow,  eta = -(f/g)*U*y
%  A + B a*a = C ln a
%  2Ba = C/a
%  and delsq eta = 2 a*a C = d/lambsq
%   so   C = d*a*a/2 lambsq 
%        B = d/(4 lambsq)
%        A = (d*a*a/lambsq)(1/2 ln a - 1/4)
%        eta = -fUg/y + 1/4 d/lambsq (r^2-a^2) + (1/2 d/lambsq) log(a)     r<a
 %            = -fUg/y + 1/2 d (a*a/lambsq)  log(r)                      r>a
 
L=0.5e6;ngrid=50;    % 1000 km box, 100 km-radius mt, 100m tall, 100 gridpoints (=2x50)
g=10; H=1000; f=1.e-4; a=1.e5; d=-100;
U=.33;
lambsq=g*H/(f*f);
Ro=U/(f*a)
Ro*H/d   %  this parameter decides whether the mountain significantly alters the flow 
%          and closed contours of psi or p or eta appear when d/H > Ro
%          roughly.
C = d*a*a/(2*lambsq);
B = d/(4*lambsq);
A = (d*a*a/lambsq)*(log(a)/2-1/4);
A
B 
C
lambsq
[x y]=meshgrid(-L:L/ngrid:L,-L:L/ngrid:L);
r=sqrt(x.*x+y.*y);
theta=atan2(y,x);   % this works in 4 quadrants whereas atan does not.
eta=-(f*U/g)*y + A + B*r.*r;
jj=find(r>a);
eta(jj)=-(f*U/g)*y(jj)+ C*log(r(jj));
contour(x,y,eta,25)
figure
meshc(x,y,eta)
title('{\eta},{\psi},p for rotating flow over a cylindrical mountain, single layer, Ro/d = 0.03')