function [xt, yt, nevals] = brownian_dynamics_with_switching(x0, y0, params, D, dt, nsave, niter, Pswitch, dx_switch, Tswitch)
% 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
%  Pswitch - probability of trying a switching move
%  dx_switch - magnitude of switching move
%  Tswitch - switching time
%
% 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;

% Keep track of average acceptance probability.
Paccept_average = 0.0;
nswitchtrials = 0.0;

% Integrate
for n = 1:niter
  if (mod(n,1000) == 0)
    n
  end
  
  % Select move type.
  if (rand > Pswitch) 
    % Brownian dynamics  

    % 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
  else
    % Switching
    
    % Generate random variables.
    Ry = randn(Tswitch,1)*sqrt(2*D*dt);   %%% random variables, variance 2D(dt)    

    % Choose sign randomly.
    sign = +1.0;
    if (rand < 0.5)
      sign = -1.0;
    end
    % DEBUG
    %sign = 2*rand-1;
    
    % x is driven
    dx = sign * dx_switch / Tswitch;
    
    % Integrate Tswitch steps, accumulating work
    xold = x0; yold = y0;
    work = 0.0;
    [potential, gradient] = grad_V(x0, y0, params);    

    if (rand < 0.5)
      % Perturb, then propagate.
      for t = 1:Tswitch
	%% Drve x.

	% Increment x.
	x0 = x0 + dx;
	% Evaluate potential and gradient, incrementing work
	work = work - potential;
	[potential, gradient] = grad_V(x0, y0, params);
	nevals = nevals + 1;
	work = work + potential;	

	%% Propagate y.

	% Compute displacement in y
	dy = Ry(t) - gradient(2)*dt;    
	% Update position in y.
	y0 = y0 + dy;
      end
    else
      % Propagate, then perturb.
      for t = 1:Tswitch
	%% Propagate y.

	% Compute displacement in y
	dy = Ry(t) - gradient(2)*dt;    
	% Update position in y.
	y0 = y0 + dy;

	%% Drive x.

	% Increment x.
	x0 = x0 + dx;
	% Evaluate potential and gradient, incrementing work
	work = work - potential;
	[potential, gradient] = grad_V(x0, y0, params);
	nevals = nevals + 1;
	work = work + potential;	
      end
    end
    
    % Accept or reject.
    Paccept = min(1.0, exp(-work));
    Paccept_average = Paccept_average + Paccept;
    nswitchtrials = nswitchtrials + 1;
    % disp(sprintf('work = %8f, Paccept = %8f', work, Paccept));
    if (rand < Paccept)
      %disp('accepted')
    else
      % restore old positions
      x0 = xold; y0 = yold;
      %disp('rejected');
    end
      
  end 
    
  % Store
  xt(n+1) = x0;
  yt(n+1) = y0;
end

Paccept_average = Paccept_average / nswitchtrials;
disp(sprintf('nonequilibrium switching trials: average acceptance %8f in %d trials', Paccept_average, nswitchtrials));

return
