% Make a plot of 2D model potential.

% Import model potential parameters.
landscapes; % load all landscapes
params = landscape.AB; % select landscape

% Generate a 2D grid to plot landscape.
xmin = -20; xmax = +20;
ymin = -20; ymax = +20;
npoints = 100;
xpoints = linspace(xmin, xmax, npoints);
ypoints = linspace(ymin, ymax, npoints);
Uxy = zeros(npoints, npoints);
for ix = 1:npoints
  for iy = 1:npoints
    x = xpoints(ix);
    y = ypoints(iy);
    [potential, gradient] = grad_V(x,y,params);
    Uxy(ix,iy) = potential;
  end
end

% Make plot
clf;

% Make contour plots in left and right panels.
contour_max = 6.0; % in kT
contour_step = 0.5; % in kT
contours = 0:contour_step:contour_max;
plot_handles = []; % storage for contour plot axis handles

ax = subplot(1,2,1);
contour(xpoints, ypoints, Uxy, contours);
axis equal;
axis([xmin xmax ymin ymax]);
colorbar;
hold on;
plot_handles = [plot_handles ax];
title('Brownian dynamics');

ax = subplot(1,2,2);
contour(xpoints, ypoints, Uxy, contours);
axis equal;
axis([xmin xmax ymin ymax]);
colorbar;
hold on;
plot_handles = [plot_handles ax];
title('nonequilibrium moves');
xlabel('driven');
ylabel('slaved');

% Generate a Brownian dynamics trajectory.
x0 = -10; y0 = -10; % initial point
D = 0.1; % diffusion constant
dt = 1; % timestep
nsave = 10; % number of steps between saving coordinates
niter = 10000; % number of nsteps periods to simulate
[xt, yt, nevals] = brownian_dynamics(x0, y0, params, D, dt, nsave, niter);

% Plot trajectory.
plot(plot_handles(1), yt, xt, '-'); % NOTE: plotting is reversed for some crazy reason!

% Count number of crossings of approximate dividing surface.
m = -0.7143; % slope
b = -5; % y-intercept
dividing_surface_x = [xmin xmax];
dividing_surface_y = m * dividing_surface_x + b;
plot(dividing_surface_x, dividing_surface_y, 'r-');
s_t = (yt > m*xt+b); % state assignments
ncrossings = sum((s_t(1:end-1) & ~s_t(2:end)) + (~s_t(1:end-1) & s_t(2:end)));
fprintf('ncrossings = %d\n', ncrossings);

% Estimate diffusion constant.
dx2_t = (xt(2:end) - xt(1:end-1)).^2 + (yt(2:end) - yt(1:end-1)).^2;
Edx2 = mean(dx2_t);
dEdx2 = std(dx2_t) / sqrt(niter);
rms = sqrt(Edx2);
drms = 0.5 / rms * dEdx2;
disp(sprintf('rms = %8.3f +- %.3f', rms, drms));

% Estimate efficiency.
efficiency = rms / nevals;
defficiency = drms / nevals;
disp(sprintf('efficiency = %8.3e +- %8.3e', efficiency, defficiency));

% Try out nonequilibrium switching.
Pswitch = 0.5; % probability of trying a switching move
dx = 8; % size of switching move
Tswitch = 10; % switching time (steps)
[xt, yt, nevals] = brownian_dynamics_with_switching(x0, y0, params, D, dt, nsave, niter, Pswitch, dx, Tswitch);

% Plot 
plot(plot_handles(2), yt, xt, '-'); % NOTE: plotting is reversed for some crazy reason!

% Estimate diffusion constant.
dx2_t = (xt(2:end) - xt(1:end-1)).^2 + (yt(2:end) - yt(1:end-1)).^2;
Edx2 = mean(dx2_t);
dEdx2 = std(dx2_t) / sqrt(niter);
rms = sqrt(Edx2);
drms = 0.5 / rms * dEdx2;
disp(sprintf('rms = %8.3f +- %.3f', rms, drms));

% Estimate efficiency.
efficiency = rms / nevals;
defficiency = drms / nevals;
disp(sprintf('efficiency = %8.3e +- %8.3e', efficiency, defficiency));

% Count number of crossings of approximate dividing surface.
m = -0.7143; % slope
b = -5; % y-intercept
dividing_surface_x = [xmin xmax];
dividing_surface_y = m * dividing_surface_x + b;
plot(dividing_surface_x, dividing_surface_y, 'r-');
s_t = (yt > m*xt+b); % state assignments
ncrossings = sum((s_t(1:end-1) & ~s_t(2:end)) + (~s_t(1:end-1) & s_t(2:end)));
fprintf('ncrossings = %d\n', ncrossings);

