%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                      %%
%%           Doswell Vortex             %%
%%                                      %%
%%         Analytical Solution          %%
%%                                      %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  The motion is a pure rotation with
%  vanishing radial velocity and 
%  tangential velocity:
%
%           V = (sech r)^2 * tanh r
%
%  The initial temperature depends only
%  on the `North-south' direction:
%
%           T0 = - tanh y
%
%  The study is purely kinematic. The
%  velocity is held constant in time.
%
%  The evolving temperature field is
%  obtained by integration the simple
%  advection equation
%
%           dT/dt = 0
%
%  The analytical solution is easily found.
%  Each fluid parcel is carried round a
%  circle with constant angular velocity
%  omega depending only on the radius.
%  So the value of T(x,y) at time t is equal
%  to its value at a point (x0,y0) displaced
%  through an angle omega*t.
%  
%  [See Doswell, C. A., 1984: A kinematic
%   analysis of frontogenesis associated
%   with a nondivergent vortex. J. Atmos.
%   Sci., 41, 1242-1248.]
%
%  [See also EIgil Kaas, Tellus, 2007 or 8]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clf
colormap('jet')

rmax = 4; 
timemax = 30;

Ntime = 1*timemax;
dtime = timemax/Ntime;

Nx = 101;             Ny = 101;
Nx = 201;             Ny = 201;

xmax = 0.7*rmax;      ymax = xmax;
dx = 2*xmax/(Nx-1);   dy = dx;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Get fields on Cartesian grid.

for nx=1:Nx
  X(nx) = -xmax+(nx-1)*dx;
end
for ny=1:Ny
  Y(ny) = -ymax+(ny-1)*dy;
end
[XX,YY] = meshgrid(X,Y);
XX = XX'; YY = YY';

for ny=1:Ny
for nx=1:Nx
  x = X(nx);
  y = Y(ny);
  r = sqrt(x^2+y^2);
  theta = atan2(y,x);
   Vrot =  (sech(r))^2 * tanh(r);
   TT(nx,ny) = -tanh(y);
   UU(nx,ny) = -Vrot*sin(theta);
   VV(nx,ny) = +Vrot*cos(theta);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Plot the initial temperature
contourf(XX,YY,TT)
axis('square'); colorbar
drawnow
pause

%%  Plot the vortex flow
quiver(XX,YY,UU,VV)
axis('square')
hold on
SS = sqrt(UU.^2+VV.^2);
contour(XX,YY,SS);
hold off
drawnow
pause
surf(XX,YY,SS)
view(0,90); shading('interp')
axis('square'); colorbar
drawnow
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Compute the solution of the
%%  advection equation for each time.

for ntime=1:Ntime
  time = ntime*dtime;
 
  %  Compute the departure point for each (x,y)
  for ny=1:Ny
  for nx=1:Nx
    x = X(nx);
    y = Y(ny);
    r = sqrt(x^2+y^2);
    theta = atan2(y,x);
    omega = ((sech(r))^2*tanh(r))/(r+eps);
    omt = omega*time;
    r0 = r;
    theta0 = theta - omt;
    x0 = r0*cos(theta0);
    y0 = r0*sin(theta0);
    TT(nx,ny) = -tanh(y0*cos(omt)-x0*sin(omt));
  end
  end

  %  Plot the temperature at each time
  surf(XX,YY,TT)
  view(0,90)
  %  view(-90,60)
  axis('square')
  axis('off')
  shading('interp')
  colorbar
  heading = ['Time = ' num2str(time) ];
  title(heading)
  drawnow
  %  pause

  day = ntime*dtime;
  if(day==1|day==3|day==5|day==10|day==20) pause; end

end   %   End of time loop

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
