
if ~exist('start_pattern')
  start_pattern = 'patch';
end

if ~exist('rk4_init')
  rk4_init = 1;
end

% INITIALISE THE VELOCITY DATA

if rk4_init
  
  disp('Using preloaded u,v data...')
  
  % h = hmean+hoi_top+hoi_top_res;
  % h = hmean+h;
  % [u,v] = uv_from_h_h(lon,lat,h);
  % u = mask.*u;
  % v = mask.*v;
  % [dx,dy] = dxdy_h(lon,lat);
  
  % divide velocity by grid metrics (see help rk4 for why we do this)
  if ndims(u)==2
      disp('Only one set of u,v - assuming steady flow')
      u3d(:,:,1) = u./dx;
      u3d(:,:,2) = u./dx;
      v3d(:,:,1) = v./dy;
      v3d(:,:,2) = v./dy;
      t(1) = 0;
      t(2) = 365*86400; 
  else
      disp('Scaling u,v by dx,dy for all time values')
      for tindex=1:size(u,1)
          u3d(:,:,tindex) = squeeze(u(tindex,:,:))./dx;
          v3d(:,:,tindex) = squeeze(v(tindex,:,:))./dy;
      end
      t = 86400*dnum(:)';
  end
  
  % make the perimeter velocity zero to stop drifters exiting the domain and 
  % causing subscripting errors in the interpolation routines
  
  u3d([1 size(u3d,1)],:,:) = 0;
  v3d([1 size(v3d,1)],:,:) = 0;
  u3d(:,[1 size(u3d,2)],:) = 0;
  v3d(:,[1 size(v3d,2)],:) = 0;
  
  % the nans on land mess things up in the interpolation
  % set u,v(land) to zero and the drifters just stall (run aground)
  u3d = change(u3d,'==',NaN,0);
  v3d = change(v3d,'==',NaN,0);
  
end

% INITIALIZE THE DRIFTERS
% need 2 row vectors of longitude and latitude

disp(['Drifter start pattern selection is: ' start_pattern])
pcolorjw(lon,lat,h(1,:,:))
hold on
plot(lon_0,lat_0,'r^')

switch start_pattern
  
  case 'click'
    disp('Using ginput to set initial drifter locations')
    disp(' Click on locations.  Press RETURN to end input...')
    [Lon0,Lat0] = ginput;
    Lat0 = Lat0'; % needs to be row vector
    Lon0 = Lon0'; % needs to be row vector
    disp(['  Initialized ' int2str(length(Lat0)) ' drifters'])
    
  case 'line'
    disp('Using ginput to set line of drifter start locations')
    disp(' Click on end points of line...')
    [Lon0,Lat0] = ginput(2);
    numd = input('How many drifters? ');
    Lon0 = interp1([1 numd],Lon0',1:numd);
    Lat0 = interp1([1 numd],Lat0',1:numd);
  
  case 'patch'
    disp('Using ginput to set patch of drifter start locations')
    disp(' Click on 2 opposite corners of patch ...')
    [Lon0,Lat0] = ginput(2);
    numdx = input('How many points in x direction? ');
    numdy = input('How many points in y direction? ');
    Lon0 = interp1([1 numdx],Lon0',1:numdx);
    Lat0 = interp1([1 numdy],Lat0',1:numdy);
    [Lon0,Lat0] = meshgrid(Lon0,Lat0);
    Lon0 = Lon0(:)';
    Lat0 = Lat0(:)';
    
  otherwise
    
end
  

% INITIALIZE EVERYTHING ELSE FOR THE PARTICLE TRACKING

[X,Y] = meshgrid(1:size(lon,2),1:size(lon,1));
x0 = interp2(lon,lat,X,Lon0,Lat0);
y0 = interp2(lon,lat,Y,Lon0,Lat0);

tp = 86400* (dnum(1):0.5:dnum(end));
xp = x0;
yp = y0;

% t = range(tp)'; % JUST A HACK FOR A STEADY VELOCITY FIELD

% INTEGRATE THE PARTICLE TRAJECTORIES

for i=2:length(tp)
  disp(['Doing ' int2str(i) ' of ' int2str(length(tp)) '; tp = ' ...
      num2str(tp(i))])
  [xp(i,:),yp(i,:)] = rk4(X,Y,t',u3d,v3d,xp(i-1,:),yp(i-1,:),...
      tp(i-1),tp(i),size(X,2),size(X,1));
end

lonp = interp2(1:size(lon,2),1:size(lon,1),lon,xp,yp);
latp = interp2(1:size(lon,2),1:size(lon,1),lat,xp,yp);

