% PUCKS_ON_ICE.M
%
% Pucks_on_ice examines the effect of an impulsively-imposed body
% force upon a series of "hockey pucks" of uniform mass that are
% distributed about the surface of a frozen ocean.  The pucks are
% initially at rest. The body force, which pushes the pucks to the
% northeast, acts for 3 days and then stops.  The integration of
% the ensuing motion continues for period of time after the forcing
% has ceased.
%
% This is an "f-plane" calculation in two dimensions, in which the
% Coriolis force is set to a fixed latitude (30 degrees north latitude 
% in this example).  There are a number of "forces" in this problem, 
% in addition to the body force.  First, Rayleigh friction is weak and 
% opposed to the motion of the particles.  If the particle velocity is u, 
% then the frictional force per unit mass is expressed as -r*u.  Second,
% because the frozen ocean surface is not flat, gravity (g, directed
% downwards) accelerates the pucks down the sloping surface.  In
% the x-direction this acceleration, or force per unit mass, is given
% by -g*(dh/dx), where dh is the change in surface height over a
% distance, dx.  In this example the inertial period, 2\pi/f, is 24 hours.
% The equations used are as follows:
%
% dx/dt = u
% dy/dt = v
% du/dt - fv = -ru - g(dh/dx) + Fx
% dv/dt + fu = -rv - g(dh/dy) + Fy,
%
% where (x,y) are the eastward and northward particle position, (u,v)
% are the eastward and northward particle velocity, f is the Coriolis
% parameter, and (Fx,Fy) are the eastward and northward components of
% the external body force.
%
% The frozen sea surface has a "bump" that is Gaussian in shape, 
% centered at (x_m,y_m), with an amplitude, h_0, and a length scale
% defined by (lx,ly).  The sea surface height is given by the following:
%
% h(x,y) = h_0*exp(-((x-x_m)/lx)^2 -((y-y_m)/ly)^2),
% where h_0 = 5 cm, x_m = y_m = 400 km, and lx = ly = 30 km.
%
% This example illustrates several important phenomena:
% (1) Particles move at right angles to the surface force (Ekman layers).
% (2) Particles oscillate with an inertial period.
% (3) Particles near the sea surface bump get "trapped" and continue to
%     move long after those far away from the bump have come to rest.
%     This illustrates the "geostrophic" balance.
%
% All units used in the calculation are given in MKS, but are
% converted to km and hours for display.  In Figure 1 the trajectories
% of 30 particles are plotted with an 'x' every 24 hours. In Figure 2
% the trajectories of two particular particles are shown (A & B) relative 
% to where they started out. 
%
% Terry Joyce, July 2002, modified for low friction limit Sept. 2003

%----------------------------------------------------------------------------
% DISPLAY TEXT

clear
close all
display([' Pucks on ice']);
str2(1)  = {'Pucks on ice, Terry Joyce, July 2002:' };
str2(2)  = {'  '};
str2(3)  = {'Pucks\_on\_ice examines the effect of an impulsively-imposed body'};
str2(4)  = {'force upon a series of "hockey pucks" of uniform mass that are'};
str2(5)  = {'distributed about the surface of a frozen ocean.  The pucks are'};
str2(6)  = {'initially at rest. The body force, which pushes the pucks to the'};
str2(7)  = {'northeast, acts for 3 days and then stops.  The integration of'};
str2(8)  = {'the ensuing motion continues for period of time after the forcing'};
str2(9)  = {'has ceased.'};
str2(10) = {''};
str2(11) = {'This is an "f-plane" calculation in two dimensions, in which the'};
str2(12) = {'Coriolis force is set to a fixed latitude (30 degrees north latitude'}; 
str2(13) = {'in this example).  There are a number of "forces" in this problem,'}; 
str2(14) = {'in addition to the body force.  First, Rayleigh friction is weak and'}; 
str2(15) = {'opposed to the motion of the particles.  If the particle velocity is u,'};
str2(16) = {'then the frictional force per unit mass is expressed as -r*u.  Second,'};
str2(17) = {'because the frozen ocean surface is not flat, gravity (g, directed'};
str2(18) = {'downwards) accelerates the pucks down the sloping surface.  In'};
str2(19) = {'the x-direction this acceleration, or force per unit mass, is given'};
str2(20) = {'by -g*(dh/dx), where dh is the change in surface height over a'};
str2(21) = {'distance, dx.  In this example the inertial period, 2\pi/f, is 24 hours.'};
str2(22) = {'The equations used are as follows:'};
str2(23) = {''};
str2(24) = {'dx/dt = u'};
str2(25) = {'dy/dt = v'};
str2(26) = {'du/dt - fv = -ru - g(dh/dx) + Fx'};
str2(27) = {'dv/dt + fu = -rv - g(dh/dy) + Fy,'};
str2(28) = {''};
str2(29) = {'where (x,y) are the eastward and northward particle position, (u,v)'};
str2(30) = {'are the eastward and northward particle velocity, f is the Coriolis'};
str2(31) = {'parameter, and (Fx,Fy) are the eastward and northward components of'};
str2(32) = {'the external body force.'};
str2(33) = {''};
str2(34) = {'The frozen sea surface has a "bump" that is Gaussian in shape, '};
str2(35) = {'centered at (x_m,y_m), with an amplitude, h_0, and a length scale'};
str2(36) = {'defined by (lx,ly).  The sea surface height is given by the following:'};
str2(37) = {''};
str2(38) = {'h(x,y) = h_0*exp(-((x-x_m)/lx)^2 -((y-y_m)/ly)^2),'};
str2(39) = {'where h_0 = 5 cm, x_m = y_m = 400 km, and lx = ly = 30 km.'};
str2(40) = {''};
str2(41) = {'This example illustrates several important phenomena:'};
str2(42) = {'(1) Particles move at right angles to the surface force (Ekman layers).'};
str2(43) = {'(2) Particles oscillate with an inertial period.'};
str2(44) = {'(3) Particles near the sea surface bump get "trapped" and continue to'};
str2(45) = {'    move long after those far away from the bump have come to rest.'};
str2(46) = {'    This illustrates the "geostrophic" balance.'};
str2(47) = {''};
str2(48) = {'*** Hit any key to begin...& be patient! ***'};

hf3 = figure(5);
clf
set(hf3,'Position',[40 40 700 900])
set(gca,'Visible','off')
text(-0.1, 0.50, str2,'FontSize',10,'HorizontalAlignment','left')
pause

%----------------------------------------------------------------------------
% START
%----------------------------------------------------------------------------
% Define parameters

count=0;				% number of particles
t1=0:1000; 			% time in hours 
t=3600*t1;			% time in seconds
rfact=.05;			% set the coefficient for Rayleigh friction
%-----------------> here are some key paprameters that can be changed
h0 = .05;			% set topography here
lat = 30;			% set latitude here
%----------------------------------------------------------------------------
% Set initial positions (x0,y0) in [m] and velocities (u0,v0) in [m/s] of 30 
% particles, then integrate the equations governing their motion (using ode45), 
% and plot each particle trajectory

more off
figure(1), clf
for i=1:5
 for j=1:6
   count =count+1;
   display([' calculating trajectory for particle # ' num2str(count),' of 30']); 
   x0(i,j)=(300+50*(i-1))*10^3;
   u0(i,j)=0;
   y0(i,j)=(300+40*(j-1))*10^3;
   v0(i,j)=0;

% The function 'ode45' integrates the differential equations in 'deriv', given the 
% integration time, t, and the initial position (x0,y0) and velocity (u0,v0) of 
% the particle, and the topography amplitude, h0 and the reference latitude.
% The last two don't change, but putting them here passes them to the function.
[T,XY]=ode45('deriv',t,[x0(i,j),y0(i,j),u0(i,j),v0(i,j),h0,lat]);

% Plot particle trajectory over time t
   plot(XY(:,1)/10^3,XY(:,2)/10^3);
   hold on
   h=plot(XY(1,1)/10^3,XY(1,2)/10^3,'ro');
   set(h,'linewidth',2)
   figure(1)
   plot(XY(1:24:length(T),1)/10^3,XY(1:24:length(T),2)/10^3,'bx')
   axis([250 550 250 550])
   axis square

% Save trajectories of 2 particles for plotting in Figure 2
   if(i==1 & j==1)
     XYA=XY;
   end
   if(i==3 & j==3)
     XYB=XY;
   end

 end
end

% Put labels on the plot
xlabel('East distance  [km]')
ylabel('North distance  [km]')
text(x0(1,1)/10^3+5,y0(1,1)/10^3+5,'A')
text(x0(3,3)/10^3+5,y0(3,3)/10^3+5,'B')

% plot force vector
quiver(260,510,35,35,0,'k')
text(255,502,'F')
title('Pucks on ice, Fig. 1', 'fontsize', 14)

%----------------------------------------------------------------------------
% Plot close-up of particles A and B in Figure 2

figure(2), clf
axis([-70 70 -70 70])
xlabel('East distance from initial position  [km]')
ylabel('North distance from initial position  [km]')
axis square
hold on
title('Pucks on ice, Fig. 2', 'fontsize', 14)
for j=1:floor(length(T)/24)
  plot((XYA(24*j-23,1)-XYA(1,1))/10^3,(XYA(24*j-23,2)-XYA(1,2))/10^3,'bx')
  plot((XYB(24*j-23,1)-XYB(1,1))/10^3,(XYB(24*j-23,2)-XYB(1,2))/10^3,'rx')
  plot((XYA(24*j-23:24*j,1)-XYA(1,1))/10^3,(XYA(24*j-23:24*j,2)-XYA(1,2))/10^3,'b-');
  plot((XYB(24*j-23:24*j,1)-XYB(1,1))/10^3,(XYB(24*j-23:24*j,2)-XYB(1,2))/10^3,'r-');
  pause(0.2)
end
legend('particle A','particle B',4)
quiver(0,0,20,20,0,'k')
text(22,22,'F')

%----------------------------------------------------------------------------
% Define and plot sea surface height/topography in the domain
% Define the topography of the sea surface to be a symmetric Gaussian "bump"
% All quantities have units of [m]

xtop=min(x0(:,1)-10^5):5*10^3:max(x0(:,1)+10^5);
ytop=xtop;
[xx,yy]=meshgrid(xtop,ytop);	% define x,y points in the domain

lx=30*10^3;                     % x length scale of bump
ly=30*10^3;                     % y length scale of bump
xm=400*10^3;                    % x location of the maximum height of bump
ym=400*10^3;                    % y location of the maximum height of bump
bump=h0*exp(-((xx-xm)/lx).^2 - ((yy-ym)/ly).^2);
hpuck=h0*exp(-((x0-xm)/lx).^2 - ((y0-ym)/ly).^2);

figure(3), clf
colormap(cool)
surf(xx/10^3,yy/10^3,bump)
shading interp
hold on
if(h0>0)
   zmin=0; zmax=2*h0;
else
   zmin=2*h0; zmax=0;
end
  
plot3(x0/10^3,y0/10^3,hpuck,'k*')
set(gca,'fontsize',14)
axis([200 600 200 600 zmin-abs(h0/10) zmax+abs(h0/10)])
view(-40,20);
%label particles A, B again
za=h0*exp(-((x0(1,1)-xm)/lx).^2 - ((y0(1,1)-ym)/ly).^2)+abs(h0/10);
zb=h0*exp(-((x0(3,3)-xm)/lx).^2 - ((y0(3,3)-ym)/ly).^2)+h0/10;
text(x0(1,1)/10^3+10,y0(1,1)/10^3+5,za,'A','fontsize',14)
text(x0(3,3)/10^3-10,y0(3,3)/10^3-5,zb,'B','fontsize',14)
%other stuff
xlabel('East distance  [km]')
ylabel('North distance  [km]')
zlabel('Height  [m]')
title('Pucks on ice, Fig. 3')

%----------------------------------------------------------------------------
% Text explaining figures

str1(1)  = {'Pucks\_on\_ice, discussion & homework problem:' };
str1(2)  = {''};
str1(3)  = {'Figure 1 shows the trajectories of 30 particles that were initially'};
str1(4)  = {'at rest. The "wind" body force, F, is applied for 3 days in the '};
str1(5)  = {'direction indicated by the vector.  Red circles indicate the initial'};
str1(6)  = {'position of each particle, and x''s mark the particle position after each'};
str1(7)  = {'day.  Many of the particles stop their translation at right angles to the'};
str1(8)  = {'wind after the forcing ceases, and their motion decays in inertial circles'};
str1(9)  = {'at a rate given by the friction.  Some of the particles continue to move'};
str1(10) = {'in a clockwise fashion around the topographic bump located at'};
str1(11) = {'(400,400) km in the plot. '};
str1(12) = {''};
str1(13) = {'Use the "zoom" feature to observe the trajectories of different particles'};
str1(14) = {'both during the forcing and after the forcing has ceased.'};
str1(15) = {''};
str1(16) = {'The trajectories of the particles labelled "A" and "B" are plotted in Figure 2.'}; 
str1(17) = {'To see these plotted again, run draw\_fig2.m.'};
str1(18) = {'Compare the difference in the motion of the two particles relative to their'};
str1(19) = {'initial positions.  Particle A is far from the topography and does not feel any'};
str1(20) = {'"force" down the slope.  Particle B, however, feels the slope, begins to'};
str1(21) = {'slide down it and is deflected to the right by the Coriolis force.  It'};
str1(22) = {'continues to move long after the forcing stops because the sloping surface'};
str1(23) = {'remains.  Except for friction (which causes the outward spiral), it is in'};
str1(24) = {'"geostrophic balance".'};
str1(25) = {''};
str1(26) = {'Figure 3 illustrates the shape of the frozen sea surface, and the initial'};
str1(27) = {'positions of the particles.'};
str1(28) = {''};
str1(29) = {'Adjust the following parameters and observe how the particle motion changes:'};
str1(30) = {'(1) The latitude, to observe how the period of the inertial oscillations changes'};
str1(31) = {'(2) The sign of h0 in "pucks\_on\_ice.m", to change the "bump" into a "dimple"'};
str1(32) = {'(3) The hemisphere, by choosing a negative latitude in "pucks\_on\_ice.m"'};
str1(33) = {''};
str1(34) = {'For each case explain what happens & show plots of the results.'};
str1(35) = {''};
str1(36) = {'Terry Joyce, 12.808 study problem, July 2002'};
hf3=figure(6);
clf
set(hf3,'Position',[40 40 700 700])
set(gca,'Visible','off')
text(-0.1, 0.50, str1,'FontSize',10,'HorizontalAlignment','left')