function [xt, yt, nevals] = brownian_dynamics(x0, y0, params, D, dt, nsave, niter)
% Generate a Brownian dynamics trajectory in a 2D model potential.
%
% ARGUMENTS
%  x0, y0 - initial coordinates
%  params - parameters for model potential
%  D - diffusion constant
%  dt - timestep
%  nsave - number of steps between saving
%  niter - number of nsave periods to simulate
%
% RETURNS
%   xt, yt - trajectory, where xt(1) is initial coordinate x0 and xt(t+1) is after t*nsave steps
%   nevals - number of potential+gradient evaluations

% Allocate storage for coordinates.
xt = zeros(nsave+1,1);
yt = zeros(nsave+1,1);

% Store initial conditions
xt(1) = x0; yt(1) = y0;

% Keep track of number of potential function evaluations.
nevals = 0;

% Integrate
for n = 1:niter
  % Generate random variables for 'nsave' steps
  Rx = randn(nsave,1)*sqrt(2*D*dt);   %%% random variables, variance 2D(dt)
  Ry = randn(nsave,1)*sqrt(2*D*dt);   %%% random variables, variance 2D(dt)  

  % Integrate nsave steps
  for t = 1:nsave
    % Evaluate gradient
    [potential, gradient] = grad_V(x0, y0, params);
    nevals = nevals + 1;
    % Compute displacement
    dx = Rx(t) - gradient(1)*dt;
    dy = Ry(t) - gradient(2)*dt;    
    % Update positions
    x0 = x0 + dx;
    y0 = y0 + dy;
  end
    
  % Store
  xt(n+1) = x0;
  yt(n+1) = y0;
end

return
