% partslope; geostrophic adjustment of a particle on a slope
%
% compute the motion of a single particle on a 
% rotating earth, numerically. This is useful for seeing how an 
% Ekman balance or a geostrophic balance appears as the time average 
% of a start up problem, or, as the instantaneous balance if there 
% is some friction and you can wait a bit. May also be useful as a 
% first look at some forms of numerical time integration of ODEs.
%
% By Jim Price, October, 1995; modified in June 2000
% Public domain for educational purposes.
%
% The equations of motion for a particle on a rotating earth are:
%      du/dt - f*v = Fx/rho - k*u,    for the east component, and
%      dv/dt + f*u = Fy/rho - k*v,    for the north component.
% 
% d()/dt is an ordinary time derivative since this model is meant
%   to describe the motion of a single particle (and not a fluid!), 
% u and v are the east and north velocity components of the particle,
% Fx and Fy are the east and north forces per unit volume applied
%   to the particle (in a fluid problem these could be boundary stress),
% rho is the nominal density of sea water,
% f is the usual Coriolis parameter, and
% k is a linear decay coefficient that represents, say, bottom drag.

clear all
close all

%  display the text string str2
str2(1) =  {'partslope:' };
str2(2) =  {'' };
str2(3) =  {'Solve for the motion of a particle sitting on a '};
str2(4) =  {'sloping sea floor and subject to Coriolis acceleration'};
str2(5) =  {'and bottom drag.  This demonstrates a local version'}; 
str2(6) =  {'of geostrophic adjustment, and the effect of '}; 
str2(7) =  {'friction on a more or less geostrophic motion. '}; 
str2(8) =  {'You may change the latitude (deg; 1 - 90), '}; 
str2(9) =  {'the e-fold time for linear drag (days; 0 - 99999), '};
str2(10) =  {'and the initial alongslope speed of the particle:'};
str2(11) =  {'   ufrac = 0 for at rest,  '};
str2(12) =  {'           = 1 for geostrophic balance, and'};
str2(13) =  {'           > 1 for super geostrophic balance).'};

hf3 = figure(9);
clf
set(hf3,'Position',[50 200 500 500])
set(gca,'Visible','off')
text(-.05, 0.70, str2,'FontSize',12,'HorizontalAlignment','left')


latitude = 20.
hlat = uicontrol('Style', 'Edit', 'Units', 'Normalized', ...
	'Position',  [0.05 0.18 0.35 0.06], ...
	'FontSize', 10, ...
	'Callback', 'eval(get(hlat,''string''))', ...
	'String', 'latitude  = 20 ');

efold = 10.;
hefold = uicontrol('Style', 'Edit', 'Units', 'Normalized', ...
	'Position',  [0.05 0.10 0.35 0.06], ...
	'Callback', 'eval(get(hefold,''string''))', ...
	'FontSize', 10, ...
	'String', 'efold = 10 ');

ufrac = 0.;
hufrac = uicontrol('Style', 'Edit', 'Units', 'Normalized', ...
	'Position',  [0.05 0.02 0.35 0.06], ...
	'Callback', 'eval(get(hufrac,''string''))', ...
	'FontSize', 10, ...
	'String', 'ufrac = 0');

hgo = uicontrol('Style', 'Pushbutton', 'Units', 'Normalized', ...
	'Position',  [0.80 0.10 0.1 0.06], ...
	'Callback', 'uiresume(9)', ...
	'FontSize', 10, ...
	'String', 'Start');

hquit = uicontrol('Style', 'Pushbutton', 'Units', 'Normalized', ...
	'Position',  [0.80 0.02 0.1 0.06], ...
	'Callback', 'lquit = 1; uiresume(9)', ...
	'FontSize', 10, ...
	'String', 'Quit');
lquit = 0;

while (lquit == 0)
disp  ('  ... waiting for control panel action.')
uiwait(9)
if lquit == 1; break; end;

day = 8.64e4;

% lat = input('Enter the latitude (deg)');
f = 2.*7.292e-5*sin(latitude*pi/180.);

if abs(f) <= 1.e-10; f = 1.e-8; end   % don't let f = 0

% The force (Fx, Fy) is evaluated as if it were a buoyanct force. 
% As now set up, the force is northward only i.e., forcex = 0 and 
% forcey > 0.

h = 100.;            % a layer thickness, m.
slope = 0.01;
g = 10.;
delrho = 0.5;
rho = 1000.;         % nominal sea water density. 

% There is a linear drag or frictional force, that will cause an 
% e-folding of an unforced current in time efold (days).

% efold = input('Enter the e-folding time for drag (days, 1 - 9999999)' );

k = 1./(efold*day);           % compute the linear drag consistent with this efold.

spd = 100.;                   % time steps per day
ndays = 10.;                  % the number of days to run (10 is plenty).
dt = (2*pi/f)/spd;            % the time step, sec. 
nstep = ndays*(2*pi/f)/dt;    % the number of time steps needed to go ndays.

u = 1:nstep; v=1:nstep; t=1:nstep;  % set up some vectors to store the solution. 

% the (geostrophic) scales used to plot things
usc = g*(delrho/rho)*slope/f;
xsc = 0.001*usc/f;
tsc = 2*pi/f;
zsc = 1000.*xsc*slope;
esc = h*usc^2;

u(1) = ufrac*usc; v(1) = 0; t(1) = 0.;      % define the initial condition. 
fricwork(1) = 0.;
%  the initial position of the particle
x(1) = 10*xsc;
y(1) = xsc;

% Begin a loop that will step the equations ahead in time.

for i=1:nstep-1

t(i+1) = t(i) + dt;      % increment the time.

Fx = 0.;
Fy = g*delrho*slope;

% A fourth order Runga-Kutta method; very accurate but relatively expensive.

du1 =  dt*( Fx/rho + f*v(i) - k*u(i) );   
dv1 =  dt*( Fy/rho - (f*u(i) + k*v(i)) ); 

du2 =  dt*( Fx/rho + f*(v(i)+dv1/2.)  - k*(u(i)+du1/2.) );
dv2 =  dt*( Fy/rho - (f*(u(i)+du1/2.) + k*(v(i)+dv1/2.)) ); 

du3 =  dt*( Fx/rho + f*(v(i)+dv2/2.)  - k*(u(i)+du2/2.) );
dv3 =  dt*( Fy/rho - (f*(u(i)+du2/2.) + k*(v(i)+dv2/2.)) );

du4 =  dt*( Fx/rho + f*(v(i)+dv3)  - k*(u(i)+du3) );
dv4 =  dt*( Fy/rho - (f*(u(i)+du3) + k*(v(i)+dv3)) );

u(i+1) = u(i) + du1/6. + du2/3. + du3/3. + du4/6.;
v(i+1) = v(i) + dv1/6. + dv2/3. + dv3/3. + dv4/6.;

forx(1,i+1) = Fx/rho;
forx(2,i+1) = f*v(i);
forx(3,i+1) = -k*u(i);

fory(1,i+1) = Fy/rho;
fory(2,i+1) =  -f*u(i);
fory(3,i+1) = -k*v(i);

x(i+1) = x(i) + dt*u(i+1)/1000.;
y(i+1) = y(i) + dt*v(i+1)/1000.;
z(i+1) = y(i)*1000*slope;

ke(i+1) = 0.5*h*(u(i+1)^2 + v(i+1)^2);
pe(i+1) = -g*delrho*h*(y(i+1)*slope);
fricwork(i+1) = fricwork(i) - k*h*(u(i+1)^2 + v(i+1)^2)*dt;
esum(i+1) = ke(i+1) +  pe(i+1) - fricwork(i+1);

end                  % the end of the time stepping loop.

% set default graphics things
set(0,'DefaultTextFontSize',12)
set(0,'DefaultAxesFontSize',12)
set(0,'DefaultAxesLineWidth',1.2)
set(0,'DefaultLineLineWidth',1.6)

%  make a movie-like plot of position and depth

x3 = 0:4:80;
y3 = 0:1:4;


nx3 = max(size(x3));
ny3 = max(size(y3));

for i=1:nx3
for j=1:ny3
	z3(j,i) = -y3(j)*xsc*slope*1000./zsc;
end
end

notthis = 1
if notthis == 0
figure(4)
clf reset
hsurf4 = surfl (x3, y3, z3);
set (hsurf4, 'FaceColor', [.8 .9 1.]);
depth = -y*slope*1000.;
hold on
hp9 = plot3(x/xsc, y/xsc, (depth+20)/zsc, 'r');
set(hp9, 'Linewidth', 1.6)
view(-70, 60);
xlabel(' x/(g`\alpha/f^2)')
ylabel(' y/(g`\alpha/f^2)')
zlabel(' z/(g`\alpha^2/f^2)')
axis([0 60 0 4 -4 0])
title('dense particle on a rotating slope w/ friction')
end


m = 0;
for j=1:20:nstep
  m = m + 1;

  if m == 1
   figure(1)
    clf reset
    hsurf1 = surfl(x3, y3, z3);
    set (hsurf1, 'FaceColor', [.8 .9 1.]);
    depth = -y*slope*1000.;
    hold on
    hp9 = plot3 (x(1)/xsc, y(1)/xsc, (depth(1)+10)/zsc, 'r*');
    set(hp9, 'erasemode', 'none')
    % Do not give initial line plot.
    % plot3 (x/xsc, y/xsc, (depth+50)/zsc, 'm-');

    view(-70, 60);
    xlabel(' x/(g`\alpha/f^2)')
    ylabel(' y/(g`\alpha/f^2)')
    zlabel(' z/(g`\alpha^2/f^2)')
    axis([0 60 0 4 -4 0])
    title('dense particle on a slope w/ rotation and friction')
    
tip = num2str(t(j)/(2*pi/f));
tipstr = [tip, ' IP'];
hth = text(50., 3.5, 0., tipstr);

    drawnow
    end

 if m > 1
    hold on
    plot3 (x(jl:j)/xsc, y(jl:j)/xsc, (depth(jl:j)+10)/zsc, 'r-', ...
          'erasemode', 'none', 'LineWidth', 1.6);
    pause(0.1)   %  slow things a bit
    if mod(m,2) == 1
      tip = num2str(t(j)/(2*pi/f));
      tipstr = ['t/(2\pi/f) = ', tip];
      set(hth, 'String', tipstr)
    end
    drawnow
 end   %  if to plot or not

  jl = j;
  end  %  loop on j for plotting each time step

figure(2)
clf reset

subplot(2,1,1)
plot(t/tsc,u/usc,t/tsc,v/usc)
grid
legend('along','across')
legend boxoff
xlabel('time/(2\pi/f)')
ylabel('U/(g`\alpha/f)')
title('along and across slope current components')

subplot(2,1,2)
plot(x/xsc, y/xsc)
grid
axis('square')
xlabel('along slope distance, x/(g`\alpha/f^2)')
ylabel('across slope, y/(g`\alpha/f^2)')
axis('equal')

figure(3)
clf reset
subplot(2,1,1)

pe1 = pe(2);
pe = pe - pe1;
esum1 = esum(2);
esum = esum - esum1;

plot(t/tsc, ke/esc, t/tsc, pe/esc, t/tsc, fricwork/esc, t/tsc, esum/esc)
xlabel('time/(2\pi/f)')
ylabel('energy/(g`\alpha/f)^2')
legend('kinetic E', 'potential E', 'frictional W', 'E + W')
legend boxoff
title('work and energy changes')
hold on
plot(t/tsc, 0*t, 'linewidth', 0.5)

format compact

% Check the transport and a few other things.
% Assumes that the external force is in the north direction only.

eky = -Fx/(rho*f)        %  Ekman numbers 
ekx = Fy/(rho*f)

ekman = k/f              %  Ekman number
alpha = atan(ekman)      %  
ufric = Fy/( rho*(f*cos(alpha) + k*sin(alpha)) )  
urat = ufric/ekx
vfric = ufric*sin(alpha)

% the end used to stop the script with the quit button 
end  % while on the time