clear all; close all; clc;
% Solves semiclassical eqns of motion for multidimensional tunnelling in W_E.
global E info keep

L = 50; W = 50; ff = 0.06; % L,W are system dimensions. For archaic reasons (that you don't need to worry about), they are in units of HALF \ell, the unit of distance used in the paper. So set them to TWICE the value you want in units of \ell. ff is the fill factor, called 'f' in the paper. The definition of this parameter is the same as in the paper, so it is easy to set.

count = 1; % Noise realisation label
folder = strcat('L',int2str(L),'W',int2str(W),'ff',num2str(ff),'P');
folder = strrep(folder,'.','_');
cd SimCity; % pick folder
cd(folder);
load(strcat('wrkspc',int2str(count),'.mat'));
cd ../..

% The above code loads a pre-saved workspace generated by ProcessNetwork.m

keep = [];
for k=1:size(info,1)
    keep = [keep; info{k,5}];
end

E = 0.3; % Set the energy of interest, in units of E_0

% Initial point for evolution can be a maximum (of u):
%x0 = maxima(7,1); y0 = maxima(7,2);

% or an arbitrary point (you need this when crossing turning surfaces).
x0 = 17.74; y0 = 16.675;

theta = 3.2; % Angle for initial velocity. THIS NEEDS GUESSING (see HowTo_Table1.txt for instructions) 
vecy = [cos(theta) sin(theta)];
e = 0.1; % Step size along the trajectory for plotting, *not* solving
[s,r] = ode45(@ClassicalEqns,(0:e:15),[x0; y0; vecy(1); vecy(2)]); % Final point along the arc length needs adjustment (at the moment, it's arbitrarily set to 15). If you run this and get an error message, reduce this number as it'll mean that your trajectory has gone so far that it's exiting the system.

figure; grid on; hold on;
% Plot the processed valley network
plot(maxima(:,1),maxima(:,2),'bx');
plot(minima(:,1),minima(:,2),'rx');
plot(saddles(:,1),saddles(:,2),'gx');

for k=1:length(C)
    traj = C{k,1};
    plot3(traj(1,:),traj(2,:),traj(3,:),'b.');
    traj = C{k,2};
    plot3(traj(1,:),traj(2,:),traj(3,:),'r.');
end
% and the candidate minimal paths from LLT
for k=1:length(C2)
    traj = C2{k,1};
    plot3(traj(1,:),traj(2,:),traj(3,:),'g.');
    traj = C2{k,2};
    plot3(traj(1,:),traj(2,:),traj(3,:),'k.');
end

% Then plot your trajectory
plot(r(:,1),r(:,2),'md');

ue = zeros(size(r,1),1);
for k=1:size(r,1)
    part = find(keep(:,1)<=r(k,1) & keep(:,2)>=r(k,1) & keep(:,3)<=r(k,2) & keep(:,4)>=r(k,2));
    part(2:end) =  [];
    u = info{part,2};
    ue(k) = feval(u,r(k,1),r(k,2));
end

% and integrate along it to get the cost according to the Agmon measure
decay = e*trapz(real(sqrt(1./ue-E)))

% This bit plots the turning surfaces for the classical trajectory -- helps a LOT!
if 1
%figure; grid on; hold on;
for part=1:size(info,1)
    
u = info{part,2};
    
if keep(part,1)==0 && keep(part,2)~=Lf
    x = (keep(part,1)+0.5:0.1:keep(part,2));
elseif keep(part,2)==Lf && keep(part,1)~=0
    x = (keep(part,1):0.1:keep(part,2)-0.5);
elseif keep(part,2)==Lf && keep(part,1)==0
    x = (keep(part,1)+0.5:0.1:keep(part,2)-0.5);
else
    x = (keep(part,1):0.1:keep(part,2));
end

if keep(part,3)==0 && keep(part,4)~=Wf
    y = (keep(part,3)+0.5:0.1:keep(part,4));
elseif keep(part,4)==Wf && keep(part,3)~=0
    y = (keep(part,3):0.1:keep(part,4)-0.5);
elseif keep(part,4)==Wf && keep(part,3)==0
    y = (keep(part,3)+0.5:0.1:keep(part,4)-0.5);
else
    y = (keep(part,3):0.1:keep(part,4));
end
    
[X,Y] = meshgrid(x,y);
ue = feval(u,X,Y);
%surface(x,y,1./ue); colormap('jet'); shading interp;
contour(X,Y,1./ue,[E,E]);

%ux = info{part,3}; uxe = feval(ux,X,Y);
%uy = info{part,4}; uye = feval(uy,X,Y);
%Vx = uxe;
%Vy = uye;
%quiver(X,Y,Vx,Vy);
end

end
