%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  This routine contains commands for plotting wave solutions to the
%  equations of motion for a shallow water model on an equatorial beta
%  plane.  This requires the function "hermite_phys.m", also available
%  from the course website.  Plotting requires "plot_wave_structure.m",
%  from the course website.
%
%  The structures follow Matsuno, 1966.  
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Dan Vimont, 30 January, 2007

clear



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Start by defining dimensional parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

RADUS = 6.375e6;            %  Radius of Earth
OMEGA = 2*pi/(60*60*24);    %  Angular velocity of Earth

beta = 2*OMEGA/RADUS;       %  Meridional derivative of the Coriolis
                            %  Parameter at the equator.

g = 9.8;                    %  Gravitational Acceleration
He = 400;                   %  Equivalent Depth (not used)
ce = sqrt(g*He);            %  Equivanelt wave speed
ae = sqrt(ce/beta);         %  Length scale; Matsuno

y = -5:.25:5;                %  Y-dimension (non-dimensional)
x = 2*pi*[0:.025:.999];      %  Y-dimension (non-dimensional)
[xx, yy] = meshgrid(x, y);
[ny, nx] = size(xx);

yexp = exp(-0.5*yy.^2);
xexp = exp(sqrt(-1)*xx);
dx = mean(diff(x));
dy = mean(diff(y));



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Kelvin wave:  \omega = k
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

v = zeros(size(xx));
u = real(yexp.*xexp);
phi = real(yexp.*xexp);
scl = 1./max(max(phi));

%  Plot results
figure(1);
orient tall;
wysiwyg;

%  Structure of the wave
subplot(2,1,1); cla;
  [cc1, hh1] = contour(xx, yy, scl*phi, [-1:.2:-.2], '--k');
  hold on;
    [cc2, hh2] = contour(xx, yy, scl*phi, [.2:.2:1], '-k');
    [cc3, hh3] = contour(xx, yy, scl*phi, [0 0], '-k');
    set(hh3, 'linewidth', 2);
    qq = quiver(xx, yy, scl*u, scl*v, 0.75);
  hold off;
  axis([0 2*pi -5 5]);
  set(gca, 'XTick', [0:.25:1]*2*pi, 'XTickLabel', 0:.25:1, 'YTick', -5:5);
  grid on;
  title('Kelvin Wave:  \phi, u, and v', 'fontsize', 12);
  
%  Convergence - centered difference
subplot(2,1,2); cla;
  dudx = ([u(:,2:end) u(:,1)] - [u(:,end) u(:,1:(end-1))])/(2*dx);
  [cc1, hh1] = contour(xx, yy, -dudx, [-1:.2:-.2], '--k');
  hold on;
    [cc2, hh2] = contour(xx, yy, -dudx, [.2:.2:1], '-k');
    [cc3, hh3] = contour(xx, yy, -dudx, [0 0], '-k');
    set(hh3, 'linewidth', 2);
  hold off;
  axis([0 2*pi -5 5]);
  set(gca, 'XTick', [0:.25:1]*2*pi, 'XTickLabel', 0:.25:1, 'YTick', -5:5);
  grid on;
  title('Kelvin Wave:  du/dx', 'fontsize', 12);

  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Rossby Waves:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = 1;                  %  Set this to the order you're interested in
k = -0.5;               %  Long Waves
k = -5;                 %  Short Waves
om = -k / (k^2 + 2*N+1);

Hn = hermite(y, N)*ones(1, nx);
Pn = yexp.*Hn;
Hnm1 = hermite(y, N-1)*ones(1, nx);
Pnm1 = yexp.*Hnm1;
Hnp1 = hermite(y, N+1)*ones(1, nx);
Pnp1 = yexp.*Hnp1;
u = (0.5*(om+k)*Pnp1 + N*(om-k)*Pnm1).*xexp;
v = (sqrt(-1)*(om^2 - k^2)*Pn).*xexp;
phi = (0.5*(om+k)*Pnp1 - N*(om-k)*Pnm1).*xexp;

u = real(u); v = real(v); phi = real(phi);
scl = 1./max(max(phi));

%  Plot results
figure(2);
colormap([0.3 1 0.8]'*[1 1 1]);
orient tall;
wysiwyg;

%  Structure of the wave
subplot(2,1,1); cla;
  [cc1, hh1] = contour(xx, yy, scl*phi, [-1:.2:-.2], '--k');
  hold on;
    [cc2, hh2] = contour(xx, yy, scl*phi, [.2:.2:1], '-k');
    [cc3, hh3] = contour(xx, yy, scl*phi, [0 0], '-k');
    set(hh3, 'linewidth', 2);
    qq = quiver(xx, yy, scl*u, scl*v, 0.75);
  hold off;
  axis([0 2*pi -5 5]);
  set(gca, 'XTick', [0:.25:1]*2*pi, 'XTickLabel', 0:.25:1, 'YTick', -5:5);
  grid on;
  title(['Rossby Wave, N = ' num2str(N) ':  \phi, u, and v'], 'fontsize', 12);
  
%  Vorticity Advection
subplot(2,1,2); cla;
  %  Calculate vorticity:
  dudy = ([u(2:end,:); u(end,:)] - [u(1,:); u(1:(end-1),:)])/(2*dy);
  dvdx = ([v(:,2:end) v(:,1)] - [v(:,end) v(:,1:(end-1))])/(2*dx);
  pv = (dvdx - dudy) - yy.*phi;
  f = y'*ones(1,nx);
  vdfdy = -v;
  scl = 1./max(max(pv));
  
  %  Add contours of the Coriolis parameter
  [cc1, hh1] = contour(xx, yy, f, [-1:.2:-.2]*5, '--k');
  hold on;
    [cc2, hh2] = contour(xx, yy, f, [.2:.2:1]*5, '-k');
    [cc3, hh3] = contour(xx, yy, f, [0 0], '-k');
    set(hh3, 'linewidth', 2);
  hold off;
  set([hh1; hh2; hh3], 'color', 0.6*[1 1 1]);  %  Grey
  
  %  Add contours of PV
  hold on;
    [cc1, hh1] = contour(xx, yy, scl*pv, [-1:.2:-.2], '--k');
    [cc2, hh2] = contour(xx, yy, scl*pv, [.2:.2:1], '-k');
  hold off;
  set([hh1; hh2], 'color', 0*[1 1 1], 'linewidth', 1);

  %  Add filled contours of vorticity advection
  hold on;
    surface(xx, yy, -20*ones(size(xx)), scl*vdfdy);
  hold off;
  caxis([-1 1]*1e-1);
 
  axis([0 2*pi -5 5]);
  set(gca, 'XTick', [0:.25:1]*2*pi, 'XTickLabel', 0:.25:1, 'YTick', -5: ...
	   5);
  set(gca, 'layer', 'top');
  grid on;
  title(['Rossby Wave, N = ' num2str(N) ':  PV (cont), PVA (shade)'], ...
         'fontsize', 12);