function [x_traj,y_traj,t,solver_name] = trajectories(x0,y0,time,x,y, ...
                                                      archive_time,solver, ...
                                                      abs_tol,rel_tol,path_function)
  
%*********************************************************************
%  Function to integrate a group of trajectories using a discrete
%  archive of two-dimensional velocities.
%*********************************************************************
%  Inputs
%*********************************************************************
%  x0           --> vector of particle initial X positions 
%  y0           --> vector of particle initial Y positions 
%  time         --> vector of times where particle positions are desired
%  x            --> vector of x values for the archive grid
%  y            --> vector of y values for the archive grid
%  archive_time --> vector of time values for the archive
%  solver       --> integer specifying solver type
%                   [1=Runge-Kutta, 2=Multi-step NDF]
%  abs_tol      --> absolute tolerance (see help for 'odeset')
%  rel_tol      --> relative tolerance (see help for 'odeset')
%*********************************************************************
%  Outputs
%*********************************************************************
%  x_traj       --> 2D array of particle x positions, with dimensions
%                   like x_traj(n_particles,n_time)
%
%  y_traj       --> 2D array of particle y positions, with dimensions
%                   like y_traj(n_particles,n_time)
%  
%  t            --> time values associated with the particle positions
%                   in x_traj and y_traj
%
%  solver_name  --> text string containing the name of the solver
%*********************************************************************
%  Setup required:  
%*********************************************************************
%  In addition to the supplied input parameters, the gridded archived
%  (u,v) velocity components must be stored in 3D global variables
%  named 'u' and 'v', like:
%  
%    global u v
%    
%  where u = u(ix,jy,ntime)  
%        v = v(ix,jy,ntime)  
%
%        ix    = number of x grid points
%        jy    = number of y grid points
%        ntime = number of discrete times in the archive
%*********************************************************************
%  UNITS !!!  
%********************************************************************* 
%  This function assumes that the length units describing the grid and
%  particle initial positions (input as [x,y] and [x0,y0]) and the time
%  units describing both the integration time and archive times (input as
%  time and archive_time) are consistent with those of the archived 
%  velocities (global 3D arrays u and v).
%
%  The output length and time units are the same as those used for
%  the inputs.  
%*********************************************************************
  
global X Y T ip

default('solver',1)
default('abs_tol',1e-8)
default('rel_tol',0)
default('path_function','archive_path')

abs_tol = max([abs_tol,eps]) ;
rel_tol = max([rel_tol,100*eps]) ;
abs_tol = 0.01 ;

switch solver
 case 1
  solver_routine = 'ode45' ;
  solver_name    = 'Runge-Kutta (ode45)' ;
 case 2
  solver_routine = 'ode15s' ;
  solver_name    = 'Multi-step NDF (ode15s)' ;
end

t1_string = [datestr(time(1),23),', ',datestr(time(1),15)] ;
t2_string = [datestr(time(end),23),', ',datestr(time(end),15)] ;

%---------------------------------------------------------------------  
%  Create 3D grids describing (x,y,time) values for each velocity in 
%  the archive.  These grids are passed as global variables to the 
%  path function for use by 'interpn'
%---------------------------------------------------------------------  

[X,Y,T] = ndgrid(x,y,archive_time) ;

%---------------------------------------------------------------------  
%  Store initial positions in a convenient form
%---------------------------------------------------------------------  

ntraj = length(x0(:)) ;    

ip = [x0(:) ; y0(:)] ;

%---------------------------------------------------------------------  
%  Integrate trajectories 
%---------------------------------------------------------------------  

ode_options = odeset('RelTol',rel_tol,'AbsTol',abs_tol,'Stats','on') ;

fprintf('\n --> Started integrating %i trajectories\n',ntraj)

fprintf('\n     Solver:                %s',solver_name)
fprintf('\n     Time interval:         %s to %s',t1_string,t2_string)
fprintf('\n     Time step:             %f [%d total steps]',time(2)- ...
        time(1),length(time))
fprintf('\n     Tolerances (ABS,REL):  (%8.2d, %8.2d)\n\n',abs_tol,rel_tol)

tic

[t,pos] = eval(sprintf('%s(''%s'',time,ip,ode_options)', ...
                       solver_routine,path_function)) ;
  
fprintf('\n     INTEGRATION FINISHED.  [Elapsed time = %i seconds]\n',...
             round(toc))
  
%---------------------------------------------------------------------  
%  Return trajectory positions in x and y arrays
%---------------------------------------------------------------------  

x_traj = pos(:,1:ntraj)' ;
y_traj = pos(:,ntraj+1:2*ntraj)' ;

%---------------------------------------------------------------------  
%  Correct for desired time vectors with only two entries 
%---------------------------------------------------------------------  

if length(time) == 2
    t = t([1 end]);
    x_traj = x_traj(:,[1 end]);
    y_traj = y_traj(:,[1 end]);
end