% DERIV.M
%   function d_dt=deriv(t,posvel)
% VARIABLES:
%   d_dt = a vector of rates of change of x, y, u, v for time, t
%        = [dx/dt dy/dt du/dt dv/dt]
%          NOTE:  Time, t, is in seconds but is evaluated hourly.  
%                 The 5th & 6th elements of d_dt are constants.
%      t = time in seconds
% posvel = a vector containing initial particle position (x0,y0) and 
%          velocity (u0,v0), and the topography amplitude (h0) &  
%          reference latitude (lat) (optional)
%        = [x0 y0 u0 v0 h0 lat]
%
% UNITS:
%       m = meters
%      km = kilometers
%       s = seconds
%      hr = hours
%     rad = radians
%     deg = degrees

function d_dt=deriv(t,posvel)

%----------------------------------------------------------------------------
% Define parameters of the problem

% Define amplitude & sign of topography and reference latitude
if(length(posvel)==4)
   h0=.05;	% set default bump (if not specified in input)
   lat=30;	% set default latitude (if not specified in input)
else 
   h0=posvel(5);	% user input or h0 from main program
   lat=posvel(6);	% user input or lat from main program
end

grav = 9.8;			% gravity in [m/s^2]
Omega = (2*pi)/(24*3600); 	% rotation rate of the earth in [rad/s]
f = 2*Omega*sin(lat*(pi/180)); 	% Coriolis parameter in [1/s]
rfact=0.05;
r=max(5e-6,rfact*abs(f));			% scale friction by f to have units of 1/s

%----------------------------------------------------------------------------
% Separate input vector into individual variables for clarity

x0=posvel(1);
y0=posvel(2);
u0=posvel(3);
v0=posvel(4);

%----------------------------------------------------------------------------
% Define the topography of the sea surface to be a symmetric Gaussian "bump"
% All quantities have units of [m]

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
rfact=.02;				% set default Rayleigh friction factor
% Compute the height and slope of the topography at the particle position
				
h = h0*exp(-((x0-xm)/lx)^2 -((y0-ym)/ly)^2); 
hx = -2*h*(x0-xm)/(lx^2); 	% x-derivative, or slope of topography (dh/dx)
hy = -2*h*(y0-ym)/(ly^2); 	% y-derivative, or slope of topography (dh/dy)

%----------------------------------------------------------------------------
% Define the (wind) forcing 

tfact=3600;			% converts hours to seconds
tmax=72*tfact; 			% length of time of forcing (3 days)
if(t>0 & t<tmax)
   Fx=.000005;			% units???
   Fy=Fx;
else
   Fx=0;
   Fy=0;
end

%----------------------------------------------------------------------------
% set values of differential array, d_dt

dxdt =  u0;			% initial velocity, u0
dydt =  v0;			% initial velocity, v0

% accel = Coriolis accel + frictional accel + pressure gradient term + forcing
dudt =  f*v0 - r*u0 - grav*hx + Fx;
dvdt = -f*u0 - r*v0 - grav*hy + Fy;

dt5=0; 				% don't want topography to change with time!
dt6=0;				% don't want reference latitude to change either!
% need to convert all units from [m/s] to [km/hr]
d_dt=[dxdt dydt dudt dvdt dt5 dt6]';
%dy=dy*(fact*tfact);
return