clear all; close all; clc;

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.

% We now have to load a pre-saved workspace, saved from running ValleyNetwork.m. The folders are labelled by system size and fill factor.
folder = strcat('L',int2str(L),'W',int2str(W),'ff',num2str(ff));
folder = strrep(folder,'.','_');
count = 1; % Noise realisation label, usually an integer from 1 to 20.
cd SimCity; % Choose the right folder
cd(folder);
load(strcat('wrkspc',int2str(count),'.mat'));
cd ../..

Evec = (0:0.1:2.5); % Energy vector in units of E_0
e = 0.01; % Precision parameter, makes sense to use the same value as in ValleyNetwork.m.
Lf = L/2; Wf = W/2; % These are now in \ell.

if size(maxima,1)==1
    figure; text(1,1,'SOS!','FontSize',20);
    error('Only one domain -- waste of time');
end

% search for stuck C trajectories
for n=1:length(C)
    for k=1:2
        traj = C{n,k};
        if size(traj,2)==10001
            figure; text(1,1,'SOS!','FontSize',20);
            error('C trajectory stuck')
        end
    end
end

% step (0): resolve for trajectories ending at max -- correct them

list = [];
for n=1:length(C)
    for k=1:2
        traj = C{n,k};
        dist2max = (traj(1,end)-maxima(:,1)).^2 + (traj(2,end)-maxima(:,2)).^2;
        dist2max = sqrt(dist2max);
        if ~isempty(find(dist2max<e))
			warning('traj ending at max identified for correction')
            list = [list n];
        end
    end
end

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

if ~isempty(list)
for j=1:length(list)
    Rind = list(j);
    for k=1:2
    x0 = saddles(Rind,1); y0 = saddles(Rind,2);
    
    traj = zeros(3,10001);
    traj(1,1) = x0; traj(2,1) = y0;
    for n = 2:10001   
        part = find(keep(:,1)<=traj(1,n-1) & keep(:,2)>=traj(1,n-1) & keep(:,3)<=traj(2,n-1) & keep(:,4)>=traj(2,n-1));
        part(2:end) =  [];
        u = info{part,2};
        if n==2
            traj(3,n-1) = feval(u,traj(1,n-1),traj(2,n-1));
        end
        ux = info{part,3};
        uy = info{part,4};
            
        theta = atan2(feval(uy,traj(1,n-1),traj(2,n-1)),feval(ux,traj(1,n-1),traj(2,n-1)));
        if n==2 && k==1
            theta = theta + pi/2;
        elseif n==2 && k==2
            theta = theta - pi/2;
        end
        traj(1,n) = traj(1,n-1) - e*cos(theta);
        traj(2,n) = traj(2,n-1) - e*sin(theta);
           
        traj(3,n) = feval(u,traj(1,n),traj(2,n));
            
        dist2min = (traj(1,n) - minima(:,1)).^2 + (traj(2,n) - minima(:,2)).^2;
        dist2min = sqrt(dist2min);
        indy = find(dist2min<e);
        dist2max = (traj(1,n) - maxima(:,1)).^2 + (traj(2,n) - maxima(:,2)).^2;
        dist2max = sqrt(dist2max);
        indy2 = find(dist2max<e);
        if ~isempty(indy) || traj(1,n)<e || traj(1,n)>Lf-e || traj(2,n)<e || traj(2,n)>Wf-e
            break;
        %elseif ~isempty(indy2)
        %    warning('Aaaaaah! Terminating at a maximum!')
        %    break;
        end
        if n==10001
            warning('Trajectory used up max length & did not reach extremum!')
        end
    end
    traj(:,n+1:end) = [];
    C{Rind,k} = traj;
    end
end
end        
clear list part u ux uy theta dist2min dist2max indy indy2 x0 y0 traj Rind

% step (0.1): check and clean up spurious saddle-minimum pairs

% find pairs
pairs = [];
for n=1:size(saddles,1)
    for k=1:size(minima,1)
        dist = sqrt((saddles(n,1)-minima(k,1)).^2 + (saddles(n,2)-minima(k,2)).^2);
        if dist<5*e
			warning('spurious saddle-minimum pairs found')
            pairs = [pairs; n k];
        end
    end
end

% find all trajectories that join up to the minima
if ~isempty(pairs)
joints = [];
for n=1:size(C,1)
    for k=1:2
        traj = C{n,k};
        for j=1:size(pairs,1)
            dist = sqrt((traj(1,end)-minima(pairs(j,2),1)).^2 + (traj(2,end)-minima(pairs(j,2),2)).^2); 
            if dist<5*e
                if j>size(joints,1)
                    joints(j,1:2) = [n k];
                else
                    blah = length(find(joints(j,:)>0));
                    joints(j,blah+1:blah+2) = [n k];
                end
            end
        end
    end
end

if isempty(joints)
    saddles(pairs(:,1),:) = [];
    C(pairs(:,1),:) = [];
    minima(pairs(:,2),:) = [];
    pairs = [];
end
end
if ~isempty(pairs)
% if we picked up any couples by mistake
notthese = [];
for j=1:size(pairs,1) 
    if length(find(joints(j,:)>0))==2
        notthese = [notthese j];
    end
end
pairs(notthese,:) = [];
joints(notthese,:) = [];
end

if ~isempty(pairs)
% case of min-sad-min-sad : join up all into major traj and chuck extrema
chucky = [];
if length(unique([joints(:,1);joints(:,3)]))<length([joints(:,1);joints(:,3)])
    for n=1:size(pairs,1)-1
        for k=n+1:size(pairs,1)
            if joints(n,1)==joints(k,1) || joints(n,1)==joints(k,3)
                shared = joints(n,1);
            elseif joints(n,3)==joints(k,1) || joints(n,3)==joints(k,3)
                shared = joints(n,3);
            else
                continue;
            end
            chucky = [chucky n k];
            r = find(pairs(:,1)==shared);
            if joints(r,1)==shared
                keeper = joints(r,3:4);
                trans = joints(r,1:2);
            elseif joints(r,3)==shared
                keeper = joints(r,1:2);
                trans = joints(r,3:4);
            end
            traj = C{keeper(1),keeper(2)};
            traj1 = C{trans(1),trans(2)};
            traj = [traj fliplr(traj1)];
            if r==n
                nxt = k;
            elseif r==k
                nxt = k;
            end
            
            if joints(nxt,1)==shared
                trans = joints(nxt,1:2);
                trans2 = joints(nxt,3:4);
            elseif joints(nxt,3)==shared
                trans = joints(nxt,3:4);
                trans2 = joints(nxt,1:2);
            end
            traj1 = C{trans(1),trans(2)};
            traj2 = C{trans2(1),trans2(2)};
            traj = [traj traj1 fliplr(traj2)];
            if trans2(2)==1
                trans = [trans2(1) 2];
            elseif trans2(2)==2
                trans = [trans2(1) 1];
            end
            traj1 = C{trans(1),trans(2)};
            traj = [traj traj1];
            C{keeper(1),keeper(2)}=traj;
        end
    end
end
if ~isempty(chucky)
saddles(pairs(chucky,1),:) = [];
C(pairs(chucky,1),:) = [];
minima(pairs(chucky,2),:) = [];
pairs(chucky,:) = [];
joints(chucky,:) = [];
end
end
clear chucky shared r keeper trans trans2 nxt

% case of a merging at the min
if ~isempty(pairs)
split = [];
for j=1:size(pairs,1) 
    if length(find(joints(j,:)>0))==6
        if pairs(j,1) == joints(j,1)
            sadtraj = 1;
            mergers = [3 5];
        elseif pairs(j,1) == joints(j,3)
            sadtraj = 3;
            mergers = [1 5];
        elseif pairs(j,1) == joints(j,5)
            sadtraj = 5;
            mergers = [1 3];
        end
        split = [split; j sadtraj mergers];
    end
end
if ~isempty(split)
for j=1:size(split,1)
    pairs(end+1,:) = pairs(split(j,1),:);
    joints(end+1,1:4) = joints(split(j,1),[split(j,2),split(j,2)+1,split(j,3),split(j,3)+1]);
    pairs(end+1,:) = pairs(split(j,1),:);
    joints(end+1,1:4) = joints(split(j,1),[split(j,2),split(j,2)+1,split(j,4),split(j,4)+1]);
end
pairs(split(:,1),:) = [];
joints(split(:,1),:) = [];
end
end
clear split mergers sadtraj
if ~isempty(pairs)
while ~any(joints(:,end))
    joints(:,end)=[];
end
if size(joints,2)>4
    figure; text(1,1,'SOS!','FontSize',20);
    error('new cluster fuck in town')
end
end
% case of sad-min, keep one sad, merge trajs
if ~isempty(pairs)
for j=1:size(pairs,1)    
    if pairs(j,1) == joints(j,1)
        traj = C{joints(j,3),joints(j,4)};
        traj1 = C{joints(j,1),joints(j,2)};
        traj = [traj fliplr(traj1)];
        if joints(j,2)==1
            traj1 = C{joints(j,1),2};
        elseif joints(j,2)==2
            traj1 = C{joints(j,1),1};
        end
        traj = [traj traj1];
        C{joints(j,3),joints(j,4)} = traj;
    elseif pairs(j,1) == joints(j,3)
        traj = C{joints(j,1),joints(j,2)};
        traj1 = C{joints(j,3),joints(j,4)};
        traj = [traj fliplr(traj1)];
        if joints(j,4)==1
            traj1 = C{joints(j,3),2};
        elseif joints(j,4)==2
            traj1 = C{joints(j,3),1};
        end
        traj = [traj traj1];
        C{joints(j,1),joints(j,2)} = traj;
    end    
end
    saddles(pairs(:,1),:) = [];
    C(pairs(:,1),:) = [];
    minima(pairs(:,2),:) = [];
end

clear pairs dist joints traj traj1 notthese blah

% step (0.25): resolve for trajectory pairs going from same saddle
% and running on top of each other.

same = [];
for n=1:size(saddles,1)-1
    for k=n+1:size(saddles,1)
        dist = sqrt((saddles(n,1)-saddles(k,1)).^2 + (saddles(n,2)-saddles(k,2)).^2);
        if dist<10*e
			warning('C traj pair overlap -- scheduled for sorting')
            same = [same; n k];
        end
    end
end

remove = [];
for j=1:length(C)
    traj1 = C{j,1};
    traj2 = C{j,2};
    sep = (traj1(1,end)-traj2(1,end)).^2 + (traj1(2,end)-traj2(2,end)).^2;
    sep = sqrt(sep);
    if sep<e
        [r,c] = find(same==j);
        if ~isempty(r)
            remove = [remove j];
            continue;
        end
        x0 = saddles(j,1); y0 = saddles(j,2);
    
        traj = zeros(3,10001);
        traj(1,1) = x0; traj(2,1) = y0;
        for n = 2:10001
            
            part = find(keep(:,1)<=traj(1,n-1) & keep(:,2)>=traj(1,n-1) & keep(:,3)<=traj(2,n-1) & keep(:,4)>=traj(2,n-1));
            part(2:end) =  [];
            u = info{part,2};
            if n==2
                traj(3,n-1) = feval(u,traj(1,n-1),traj(2,n-1));
            end
            ux = info{part,3};
            uy = info{part,4};
            
            theta = atan2(feval(uy,traj(1,n-1),traj(2,n-1)),feval(ux,traj(1,n-1),traj(2,n-1)));
            if n==2
                traj(1,n) = traj1(1,1) - (traj1(1,3) - traj1(1,1));
                traj(2,n) = traj1(2,1) - (traj1(2,3) - traj1(2,1));
            else
                traj(1,n) = traj(1,n-1) - e*cos(theta);
                traj(2,n) = traj(2,n-1) - e*sin(theta);
            end
            
            traj(3,n) = feval(u,traj(1,n),traj(2,n));
            
            dist2min = (traj(1,n) - minima(:,1)).^2 + (traj(2,n) - minima(:,2)).^2;
            dist2min = sqrt(dist2min);
            indy = find(dist2min<e);
            dist2max = (traj(1,n) - maxima(:,1)).^2 + (traj(2,n) - maxima(:,2)).^2;
            dist2max = sqrt(dist2max);
            indy2 = find(dist2max<e);
            if ~isempty(indy) || traj(1,n)<e || traj(1,n)>Lf-e || traj(2,n)<e || traj(2,n)>Wf-e
                break;
            elseif ~isempty(indy2)
                warning('Aaaaaah! Terminating at a maximum!')
                break;
            end
            if n==10001
                warning('Trajectory used up max length & did not reach extremum!')
            end
        end
        traj(:,n+1:end) = [];
        C{j,2} = traj;
    end
end
saddles(remove,:)=[];
C(remove,:) = [];
        
clear remove traj1 traj2 sep part u ux uy dist2min dist2max indy indy2 x0 y0 traj 
clear dist same r c

% step (0.5): filter maxima to eliminate possible repeats

same = [];
for n=1:size(maxima,1)-1
    for k=n+1:size(maxima,1)
        dist = sqrt((maxima(n,1)-maxima(k,1)).^2 + (maxima(n,2)-maxima(k,2)).^2);
        if dist<20*e
            same = [same; n k];
        end
    end
end

clear dist

% step (0.75): eliminate false double-up minima from your list.

repeats = [];
for n=1:size(minima,1)-1
    for k=n+1:size(minima,1)
        dist = sqrt((minima(n,1)-minima(k,1)).^2 + (minima(n,2)-minima(k,2)).^2);
        if dist<10*e
			warning('double min picked up')
            repeats = [repeats; n k];
        end
    end
end

if ~isempty(repeats)
belong = cell(size(repeats));
for n=1:length(C)
    for k=1:2
        traj = C{n,k};
        if ~(traj(1,end)<e || traj(1,end)>Lf-e || traj(2,end)<e || traj(2,end)>Wf-e)
            dist2min = (traj(1,end)-minima(:,1)).^2 + (traj(2,end)-minima(:,2)).^2;
            dist2min = sqrt(dist2min);
            blah = find(dist2min<e);
            if length(blah)>1
                [~,theone] = min(dist2min(blah));
                blah = blah(theone);
            end
            [r,c]=find(repeats==blah);
            if ~isempty(r)
                for j=1:length(r)
                    matty = belong{r(j),c(j)};
                    belong{r(j),c(j)} = [matty; n k];
                end
            end
        end
    end
end

repvec = zeros(size(repeats,1),1); 
for n=1:size(repeats,1)
    matty1 = belong{n,1};
    num1 = size(matty1,1);
    matty2 = belong{n,2};
    num2 = size(matty2,2);
    if num1>=num2
        xtra = [minima(repeats(n,1),1); minima(repeats(n,1),2);0];
        for k=1:size(matty2)
            traj = C{matty2(k,1),matty2(k,2)};
            traj = [traj xtra];
            C{matty2(k,1),matty2(k,2)} = traj;
        end
        repvec(n) = repeats(n,2);
    else
        xtra = [minima(repeats(n,2),1); minima(repeats(n,2),2);0];
        for k=1:size(matty1)
            traj = C{matty1(k,1),matty1(k,2)};
            traj = [traj xtra];
            C{matty1(k,1),matty1(k,2)} = traj;
        end
        repvec(n) = repeats(n,1);
    end
end
minima(repvec,:) = [];
end
clear dist blah repeats theone repvec r c matty matty1 matty2 num1 num2 xtra belong

% step (1): eliminate "stickers" from your trajectory list.

clean = 1;
while clean == 1
Terminators = zeros(size(C));
for n=1:length(C)
    for k=1:2
        traj = C{n,k};
        if traj(1,end)<e || traj(1,end)>Lf-e || traj(2,end)<e || traj(2,end)>Wf-e
            Terminators(n,k) = 0;
        else
            dist2min = (traj(1,end)-minima(:,1)).^2 + (traj(2,end)-minima(:,2)).^2;
            dist2min = sqrt(dist2min);
            Terminators(n,k) = find(dist2min<e);
         end
    end
end

howmany = zeros(length(minima),1);
pickup = [];
for n=1:length(minima)
    howmany(n) = length(find(Terminators==n));
    if howmany(n)==1
        [r,c] = find(Terminators==n);
        pickup = [pickup; r,c];
    end
end

if isempty(find(howmany==1))
    break;
end

minima(find(howmany==1),:) = [];
C(pickup(:,1),:) = [];
saddles(pickup(:,1),:) = [];
end

clear clean Terminators traj dist2min howmany pickup r c

% step (2): connect saddles to maxima & thus ID domain boundaries & domain links.

owners = zeros(size(C));

C2 = cell(length(saddles),2);
for k=1:length(saddles)
    x0 = saddles(k,1); y0 = saddles(k,2);
    for j=1:2
        traj = zeros(3,10001);
        traj(1,1) = x0; traj(2,1) = y0;
        for n = 2:10001
            
            part = find(keep(:,1)<=traj(1,n-1) & keep(:,2)>=traj(1,n-1) & keep(:,3)<=traj(2,n-1) & keep(:,4)>=traj(2,n-1));
            part(2:end) =  [];
            u = info{part,2};
            if n==2
                traj(3,n-1) = feval(u,traj(1,n-1),traj(2,n-1));
            end
            ux = info{part,3};
            uy = info{part,4};
            
            theta = atan2(feval(uy,traj(1,n-1),traj(2,n-1)),feval(ux,traj(1,n-1),traj(2,n-1)));
            if n==2
                theta = theta + pi/2;
            end
            
            if j==2 && n==2
                traj(1,n) = traj(1,n-1) - e*cos(theta);
                traj(2,n) = traj(2,n-1) - e*sin(theta);
            else
                traj(1,n) = traj(1,n-1) + e*cos(theta);
                traj(2,n) = traj(2,n-1) + e*sin(theta);
            end
            
            traj(3,n) = feval(u,traj(1,n),traj(2,n));
            
            dist2min = (traj(1,n) - minima(:,1)).^2 + (traj(2,n) - minima(:,2)).^2;
            dist2min = sqrt(dist2min);
            indy = find(dist2min<e);
            dist2max = (traj(1,n) - maxima(:,1)).^2 + (traj(2,n) - maxima(:,2)).^2;
            dist2max = sqrt(dist2max);
            indy2 = find(dist2max<e);
            if ~isempty(indy2)
                owners(k,j)= indy2(1); % if it's a doubly-picked up max & both pass
                break;
            elseif ~isempty(indy) || traj(1,n)<e || traj(1,n)>Lf-e || traj(2,n)<e || traj(2,n)>Wf-e
                warning('Aaaaaah! Terminating at a minimum or flying out to space!')
                k
                j
                break;
            end
            if n==10001
                warning('Trajectory used up max length & did not reach extremum!')
                k
                j
            end
        end
        traj(:,n+1:end) = [];
        C2{k,j} = traj;
    end
end

clear x0 y0 traj part u ux uy theta dist2min indy dist2max indy2

% step (2.25): resolve for trajectory pairs going from same saddle
% and running on top of each other.

for j=1:length(C2)
    traj1 = C2{j,1};
    traj2 = C2{j,2};
    sep = (traj1(1,end)-traj2(1,end)).^2 + (traj1(2,end)-traj2(2,end)).^2;
    sep = sqrt(sep);
    if sep<e
		warning('trying to separate two identical C2 trajectories')
        x0 = saddles(j,1); y0 = saddles(j,2);
    
        traj = zeros(3,10001);
        traj(1,1) = x0; traj(2,1) = y0;
        for n = 2:10001
            part = find(keep(:,1)<=traj(1,n-1) & keep(:,2)>=traj(1,n-1) & keep(:,3)<=traj(2,n-1) & keep(:,4)>=traj(2,n-1));
            part(2:end) =  [];
            u = info{part,2};
            if n==2
                traj(3,n-1) = feval(u,traj(1,n-1),traj(2,n-1));
            end
            ux = info{part,3};
            uy = info{part,4};
            
            theta = atan2(feval(uy,traj(1,n-1),traj(2,n-1)),feval(ux,traj(1,n-1),traj(2,n-1)));
            if n==2
                traj(1,n) = traj1(1,1) - (traj1(1,3) - traj1(1,1));
                traj(2,n) = traj1(2,1) - (traj1(2,3) - traj1(2,1));
            else
                traj(1,n) = traj(1,n-1) + e*cos(theta);
                traj(2,n) = traj(2,n-1) + e*sin(theta);
            end
            
            traj(3,n) = feval(u,traj(1,n),traj(2,n));
            
            dist2min = (traj(1,n) - minima(:,1)).^2 + (traj(2,n) - minima(:,2)).^2;
            dist2min = sqrt(dist2min);
            indy = find(dist2min<e);
            dist2max = (traj(1,n) - maxima(:,1)).^2 + (traj(2,n) - maxima(:,2)).^2;
            dist2max = sqrt(dist2max);
            indy2 = find(dist2max<e);
            if ~isempty(indy2)
                owners(j,2)= indy2(1); % if it's a doubly-picked up max & both pass
                break;
            elseif ~isempty(indy) || traj(1,n)<e || traj(1,n)>Lf-e || traj(2,n)<e || traj(2,n)>Wf-e
                warning('Aaaaaah! Terminating at a minimum or flying out to space!')
                break;
            end
            if n==10001
                warning('Trajectory used up max length & did not reach extremum!')
            end
        end
        traj(:,n+1:end) = [];
        C2{j,2} = traj;
    end
end
        
clear traj1 traj2 sep part u ux uy dist2min dist2max indy indy2 x0 y0 traj

% in case owner not picked up at all.
[rs,cs]=find(owners==0);
if ~isempty(rs)
warning('zero owner found and search radius increased')
for k=1:length(rs)
        traj = C2{rs(k),cs(k)};
        n = size(traj,2);
        dist2max = (traj(1,n) - maxima(:,1)).^2 + (traj(2,n) - maxima(:,2)).^2;
        dist2max = sqrt(dist2max);
        indy2 = find(dist2max<20*e);
        if ~isempty(indy2)
            owners(rs(k),cs(k))= indy2(1); % if it's a doubly-picked up max & both pass
        end
end
end
clear traj dist2max indy2 rs cs

% test for both owners being the same. this indicates both C2 trajs ended up in the same place.
rown = find(owners(:,1)==owners(:,2));
if ~isempty(rown)
    figure; text(1,1,'SOS!','FontSize',20);
	error('both owners of a saddle are the same')
end

% step (0.5) continued:  eliminate possible repeated maxima

remove = [];
for k=1:size(same,1)
    [r1,c1] = find(owners==same(k,1));
    [r2,c2] = find(owners==same(k,2));
    if length(r1)>=length(r2)
        for n=1:length(r2)
            owners(r2(n),c2(n)) = same(k,1);
        end
        remove = [remove; same(k,2)];
    elseif length(r1)<length(r2)
        for n=1:length(r1)
            owners(r1(n),c1(n)) = same(k,2);
        end
        remove = [remove; same(k,1)];
    %else
    %    error('Dunno which max to keep!')
    end
end

if ~isempty(same)
    retained = unique(owners);
    maxima(remove,:) = [];
    for k=1:length(retained)
        [r,c] = find(owners==retained(k));
        for n=1:length(r)
            owners(r(n),c(n)) = k;
        end
    end
end

clear remove r1 r2 c1 c2 same retained r c

% step (2.5): search for and eliminate falsely picked up duplicate saddles
% and associated trajectories.

same = [];
for n=1:size(saddles,1)-1
    for k=n+1:size(saddles,1)
        dist = sqrt((saddles(n,1)-saddles(k,1)).^2 + (saddles(n,2)-saddles(k,2)).^2);
        if dist<10*e
			warning('found doubley picked up saddle')
            same = [same; n k];
        end
    end
end

if ~isempty(same)
    remove = [];
    for n=1:size(same,1)
        
        I1 = same(n,1);
        I2 = same(n,2);
        
        % if both saddles have 2 correct trajs but same owner for each pair
        own = [owners(I1,:); owners(I2,:)];
        traj1 = C{I1,1};
        traj2 = C{I1,2};
        endpts = [traj1(1,end), traj1(2,end); traj2(1,end), traj2(2,end)];
        traj1 = C{I2,1};
        traj2 = C{I2,2};
        endpts = [endpts; traj1(1,end), traj1(2,end); traj2(1,end), traj2(2,end)];
        
        dist1 = sqrt((endpts(1,1)-endpts(3:4,1)).^2 + (endpts(1,2)-endpts(3:4,2)).^2);
        dist2 = sqrt((endpts(2,1)-endpts(3:4,1)).^2 + (endpts(2,2)-endpts(3:4,2)).^2);
        
        [val1,ind1] = min(dist1);
        [val2,ind2] = min(dist2);
        
        if val1<e && val2<e && own(1,1)==own(1,2) && own(2,1)==own(2,2)
            
            if ind2==1
                traj = traj1;
            elseif ind2==2
                traj = traj2;
            end
            
            traj(1,1) = saddles(I1,1);
            traj(2,1) = saddles(I1,2);
            C{I1,2} = traj;
            
            traj = C2{I2,ind2};
            traj(1,1) = saddles(I1,1);
            traj(2,1) = saddles(I1,2);
            C2{I1,2} = traj;
            
            owners(I1,2) = owners(I2,ind2);
            
            remove = [remove I2];
            continue;
        elseif val1<e && val2<e && own(1,1)==own(1,2) && own(2,1)~=own(2,2)      
            remove = [remove I1];
            continue;
        elseif val1<e && val2<e && own(1,1)~=own(1,2) && own(2,1)==own(2,2)      
            remove = [remove I2];
            continue;
        elseif val1<e && val2<e && own(1,1)~=own(1,2) && own(2,1)~=own(2,2)      
            remove = [remove I2];
            continue;
        end
        
        % test first saddles
        traj1 = C2{I1,1};
        traj2 = C2{I1,2};
        if size(traj1,2)<30
            blah = size(traj1,2);
            dist1 = sqrt((traj1(1,blah-10)-traj2(1,blah-10)).^2 + (traj1(2,blah-10)-traj2(2,blah-10)).^2);
            dist2 = sqrt((traj1(1,blah)-traj2(1,blah)).^2 + (traj1(2,blah)-traj2(2,blah)).^2);
        elseif size(traj2,2)<30
            blah = size(traj2,2);
            dist1 = sqrt((traj1(1,blah-10)-traj2(1,blah-10)).^2 + (traj1(2,blah-10)-traj2(2,blah-10)).^2);
            dist2 = sqrt((traj1(1,blah)-traj2(1,blah)).^2 + (traj1(2,blah)-traj2(2,blah)).^2);
        else
            dist1 = sqrt((traj1(1,20)-traj2(1,20)).^2 + (traj1(2,20)-traj2(2,20)).^2);
            dist2 = sqrt((traj1(1,30)-traj2(1,30)).^2 + (traj1(2,30)-traj2(2,30)).^2);
        end
        if abs(dist1-dist2)<e
            remove = [remove I1];
        end
        
        % test second saddle
        traj1 = C2{I2,1};
        traj2 = C2{I2,2};
        if size(traj1,2)<30
            blah = size(traj1,2);
            dist1 = sqrt((traj1(1,blah-10)-traj2(1,blah-10)).^2 + (traj1(2,blah-10)-traj2(2,blah-10)).^2);
            dist2 = sqrt((traj1(1,blah)-traj2(1,blah)).^2 + (traj1(2,blah)-traj2(2,blah)).^2);
        elseif size(traj2,2)<30
            blah = size(traj2,2);
            dist1 = sqrt((traj1(1,blah-10)-traj2(1,blah-10)).^2 + (traj1(2,blah-10)-traj2(2,blah-10)).^2);
            dist2 = sqrt((traj1(1,blah)-traj2(1,blah)).^2 + (traj1(2,blah)-traj2(2,blah)).^2);
        else
            dist1 = sqrt((traj1(1,20)-traj2(1,20)).^2 + (traj1(2,20)-traj2(2,20)).^2);
            dist2 = sqrt((traj1(1,30)-traj2(1,30)).^2 + (traj1(2,30)-traj2(2,30)).^2);
        end
        if abs(dist1-dist2)<e
            remove = [remove I2];
        end
    end
    
    saddles(remove,:) = [];
    C(remove,:) = [];
    C2(remove,:) = [];
    owners(remove,:) = [];
end

clear dist same I1 I2 traj1 traj2 dist1 dist2 remove blah
clear own endpts dist1 dist2 val1 val2 ind1 ind2

% step (1) repeated: eliminate "stickers" from your trajectory list,
% if missed on round one due to double saddle.

clean = 1;
while clean == 1
Terminators = zeros(size(C));
for n=1:length(C)
    for k=1:2
        traj = C{n,k};
        if traj(1,end)<e || traj(1,end)>Lf-e || traj(2,end)<e || traj(2,end)>Wf-e
            Terminators(n,k) = 0;
        else
            dist2min = (traj(1,end)-minima(:,1)).^2 + (traj(2,end)-minima(:,2)).^2;
            dist2min = sqrt(dist2min);
            Terminators(n,k) = find(dist2min<e);
         end
    end
end

howmany = zeros(length(minima),1);
pickup = [];
for n=1:length(minima)
    howmany(n) = length(find(Terminators==n));
    if howmany(n)==1
        [r,c] = find(Terminators==n);
        pickup = [pickup; r,c];
    end
end

if isempty(find(howmany==1))
    break;
end

minima(find(howmany==1),:) = [];
C(pickup(:,1),:) = [];
C2(pickup(:,1),:) = [];
saddles(pickup(:,1),:) = [];
owners(pickup(:,1),:) = [];
end

clear clean Terminators traj dist2min howmany pickup r c

% step (3): complete boundary domains & ID their labels.

Lexit = []; Rexit = []; Texit = []; Bexit = [];
Lloc = []; Rloc = []; Tloc = []; Bloc = [];
for n=1:length(C)
    for k=1:2
        traj = C{n,k};
        if traj(1,end)<e
            Lexit = [Lexit traj(2,end)];
            Lloc = [Lloc; n k];
        elseif traj(1,end)>Lf-e
            Rexit = [Rexit traj(2,end)];
            Rloc = [Rloc; n k];
        elseif traj(2,end)<e
            Bexit = [Bexit traj(1,end)];
            Bloc = [Bloc; n k];
        elseif traj(2,end)>Wf-e
            Texit = [Texit traj(1,end)];
            Tloc = [Tloc; n k];
        end
    end
end

[Rexit,indy] = sort(Rexit,'ascend');
Rloc = Rloc(indy,:);
suspects = find(abs(diff(Rexit))<1e-3);
cluster = [];
for n=1:length(suspects)
    if ~isempty(intersect(suspects(n),cluster))
        continue
    end
    cluster = [suspects(n)];
    for k = n+1:length(suspects)
        if suspects(k) == suspects(k-1)+1
            cluster = [cluster suspects(k)];
        else
            break;
        end
    end
    cluster = [cluster cluster(end)+1];
    if length(cluster)==2
        own1 = owners(Rloc(cluster(1),1),:);
        own2 = owners(Rloc(cluster(2),1),:);
        common = intersect(own1,own2);
        only1 = setdiff(own1,common);
        only2 = setdiff(own2,common);
        %yc = maxima(common,2);
        y1 = maxima(only1,2);
        y2 = maxima(only2,2);
        if y2<y1 %y2<yc && yc<y1
            tempy = Rexit(cluster(2));
            Rexit(cluster(2)) = Rexit(cluster(1));
            Rexit(cluster(1)) = tempy;
            tempy = Rloc(cluster(2),:);
            Rloc(cluster(2),:) = Rloc(cluster(1),:);
            Rloc(cluster(1),:) = tempy;
        end
        Rexit(cluster) = round(Rexit(cluster)*10^3)/10^3 + (0:0.001:0.001*(length(cluster)-1));
    else
        ClusterOwners = owners(Rloc(cluster,1),:);    
        ClusterMaxPos = zeros(size(ClusterOwners));
        for k=1:length(cluster)
            ClusterMaxPos(k,1) = maxima(ClusterOwners(k,1),2);
            ClusterMaxPos(k,2) = maxima(ClusterOwners(k,2),2);
        end
        try
        [blah,r] = min(ClusterMaxPos);
        [blah2,c] = min(blah);
        order = cluster(r(c));
        common = setdiff(ClusterOwners(r(c),:),ClusterOwners(r(c),c));
        ClusterOwners(r(c),:) = [0,0];
        for k=2:length(cluster)
            [r,c] = find(ClusterOwners==common);
            order = [order cluster(r)];
            common = setdiff(ClusterOwners(r,:),ClusterOwners(r,c));
            ClusterOwners(r,:) = [0,0];
        end
        catch
        ClusterOwners = owners(Rloc(cluster,1),:);    
        [blah,r] = max(ClusterMaxPos);
        [blah2,c] = max(blah);
        order = cluster(r(c));
        common = setdiff(ClusterOwners(r(c),:),ClusterOwners(r(c),c));
        ClusterOwners(r(c),:) = [0,0];
        for k=2:length(cluster)
            [r,c] = find(ClusterOwners==common);
            order = [order cluster(r)];
            common = setdiff(ClusterOwners(r,:),ClusterOwners(r,c));
            ClusterOwners(r,:) = [0,0];
        end
        order = fliplr(order);
        end
        Rexit(cluster) = Rexit(order);
        Rloc(cluster,:) = Rloc(order,:);
        Rexit(cluster) = round(Rexit(cluster)*10^3)/10^3 + (0:0.001:0.001*(length(cluster)-1));
    end
end

[Texit,indy] = sort(Texit,'descend');
Tloc = Tloc(indy,:);
suspects = find(abs(diff(Texit))<1e-3);
cluster = [];
for n=1:length(suspects)
    if ~isempty(intersect(suspects(n),cluster))
        continue
    end
    cluster = [suspects(n)];
    for k = n+1:length(suspects)
        if suspects(k) == suspects(k-1)+1
            cluster = [cluster suspects(k)];
        else
            break;
        end
    end
    cluster = [cluster cluster(end)+1];
    if length(cluster)==2
        own1 = owners(Tloc(cluster(1),1),:);
        own2 = owners(Tloc(cluster(2),1),:);
        common = intersect(own1,own2);
        only1 = setdiff(own1,common);
        only2 = setdiff(own2,common);
        %xc = maxima(common,1);
        x1 = maxima(only1,1);
        x2 = maxima(only2,1);
        if x2>x1 %x2>xc && xc>x1
            tempy = Texit(cluster(2));
            Texit(cluster(2)) = Texit(cluster(1));
            Texit(cluster(1)) = tempy;
            tempy = Tloc(cluster(2),:);
            Tloc(cluster(2),:) = Tloc(cluster(1),:);
            Tloc(cluster(1),:) = tempy;
        end
        Texit(cluster) = round(Texit(cluster)*10^3)/10^3 - (0:0.001:0.001*(length(cluster)-1));
    else    
        ClusterOwners = owners(Tloc(cluster,1),:);    
        ClusterMaxPos = zeros(size(ClusterOwners));
        for k=1:length(cluster)
            ClusterMaxPos(k,1) = maxima(ClusterOwners(k,1),1);
            ClusterMaxPos(k,2) = maxima(ClusterOwners(k,2),1);
        end
        try
        [blah,r] = max(ClusterMaxPos);
        [blah2,c] = max(blah);
        order = cluster(r(c));
        common = setdiff(ClusterOwners(r(c),:),ClusterOwners(r(c),c));
        ClusterOwners(r(c),:) = [0,0];
        for k=2:length(cluster)
            [r,c] = find(ClusterOwners==common);
            order = [order cluster(r)];
            common = setdiff(ClusterOwners(r,:),ClusterOwners(r,c));
            ClusterOwners(r,:) = [0,0];
        end
        catch
        ClusterOwners = owners(Tloc(cluster,1),:);
        [blah,r] = min(ClusterMaxPos);
        [blah2,c] = min(blah);
        order = cluster(r(c));
        common = setdiff(ClusterOwners(r(c),:),ClusterOwners(r(c),c));
        ClusterOwners(r(c),:) = [0,0];
        for k=2:length(cluster)
            [r,c] = find(ClusterOwners==common);
            order = [order cluster(r)];
            common = setdiff(ClusterOwners(r,:),ClusterOwners(r,c));
            ClusterOwners(r,:) = [0,0];
        end
        order = fliplr(order);
        end
        Texit(cluster) = Texit(order);
        Tloc(cluster,:) = Tloc(order,:);
        Texit(cluster) = round(Texit(cluster)*10^3)/10^3 - (0:0.001:0.001*(length(cluster)-1));
    end
end

[Lexit,indy] = sort(Lexit,'descend');
Lloc = Lloc(indy,:);
suspects = find(abs(diff(Lexit))<1e-3);
cluster = [];
for n=1:length(suspects)
    if ~isempty(intersect(suspects(n),cluster))
        continue
    end
    cluster = [suspects(n)];
    for k = n+1:length(suspects)
        if suspects(k) == suspects(k-1)+1
            cluster = [cluster suspects(k)];
        else
            break;
        end
    end
    cluster = [cluster cluster(end)+1];
    if length(cluster)==2
        own1 = owners(Lloc(cluster(1),1),:);
        own2 = owners(Lloc(cluster(2),1),:);
        common = intersect(own1,own2);
        only1 = setdiff(own1,common);
        only2 = setdiff(own2,common);
        %yc = maxima(common,2);
        y1 = maxima(only1,2);
        y2 = maxima(only2,2);
        if y2>y1 %y2>yc && yc>y1
            tempy = Lexit(cluster(2));
            Lexit(cluster(2)) = Lexit(cluster(1));
            Lexit(cluster(1)) = tempy;
            tempy = Lloc(cluster(2),:);
            Lloc(cluster(2),:) = Lloc(cluster(1),:);
            Lloc(cluster(1),:) = tempy;
        end
        Lexit(cluster) = round(Lexit(cluster)*10^3)/10^3 - (0:0.001:0.001*(length(cluster)-1));
    else
        ClusterOwners = owners(Lloc(cluster,1),:);    
        ClusterMaxPos = zeros(size(ClusterOwners));
        for k=1:length(cluster)
            ClusterMaxPos(k,1) = maxima(ClusterOwners(k,1),2);
            ClusterMaxPos(k,2) = maxima(ClusterOwners(k,2),2);
        end
        try
        [blah,r] = max(ClusterMaxPos);
        [blah2,c] = max(blah);
        order = cluster(r(c));
        common = setdiff(ClusterOwners(r(c),:),ClusterOwners(r(c),c));
        ClusterOwners(r(c),:) = [0,0];
        for k=2:length(cluster)
            [r,c] = find(ClusterOwners==common);
            order = [order cluster(r)];
            common = setdiff(ClusterOwners(r,:),ClusterOwners(r,c));
            ClusterOwners(r,:) = [0,0];
        end
        catch
        ClusterOwners = owners(Lloc(cluster,1),:);
        [blah,r] = min(ClusterMaxPos);
        [blah2,c] = min(blah);
        order = cluster(r(c));
        common = setdiff(ClusterOwners(r(c),:),ClusterOwners(r(c),c));
        ClusterOwners(r(c),:) = [0,0];
        for k=2:length(cluster)
            [r,c] = find(ClusterOwners==common);
            order = [order cluster(r)];
            common = setdiff(ClusterOwners(r,:),ClusterOwners(r,c));
            ClusterOwners(r,:) = [0,0];
        end
        order = fliplr(order);
        end
        Lexit(cluster) = Lexit(order);
        Lloc(cluster,:) = Lloc(order,:);
        Lexit(cluster) = round(Lexit(cluster)*10^3)/10^3 - (0:0.001:0.001*(length(cluster)-1));
    end
end

[Bexit,indy] = sort(Bexit,'ascend');
Bloc = Bloc(indy,:);
suspects = find(abs(diff(Bexit))<1e-3);
cluster = [];
for n=1:length(suspects)
    if ~isempty(intersect(suspects(n),cluster))
        continue
    end
    cluster = [suspects(n)];
    for k = n+1:length(suspects)
        if suspects(k) == suspects(k-1)+1
            cluster = [cluster suspects(k)];
        else
            break;
        end
    end
    cluster = [cluster cluster(end)+1];
    if length(cluster)==2
        own1 = owners(Bloc(cluster(1),1),:);
        own2 = owners(Bloc(cluster(2),1),:);
        common = intersect(own1,own2);
        only1 = setdiff(own1,common);
        only2 = setdiff(own2,common);
        %xc = maxima(common,1);
        x1 = maxima(only1,1);
        x2 = maxima(only2,1);
        if x2<x1 %x2<xc && xc<x1
            tempy = Bexit(cluster(2));
            Bexit(cluster(2)) = Bexit(cluster(1));
            Bexit(cluster(1)) = tempy;
            tempy = Bloc(cluster(2),:);
            Bloc(cluster(2),:) = Bloc(cluster(1),:);
            Bloc(cluster(1),:) = tempy;
        end
        Bexit(cluster) = round(Bexit(cluster)*10^3)/10^3 + (0:0.001:0.001*(length(cluster)-1));
    else
        ClusterOwners = owners(Bloc(cluster,1),:);    
        ClusterMaxPos = zeros(size(ClusterOwners));
        for k=1:length(cluster)
            ClusterMaxPos(k,1) = maxima(ClusterOwners(k,1),1);
            ClusterMaxPos(k,2) = maxima(ClusterOwners(k,2),1);
        end
        try
        [blah,r] = min(ClusterMaxPos);
        [blah2,c] = min(blah);
        order = cluster(r(c));
        common = setdiff(ClusterOwners(r(c),:),ClusterOwners(r(c),c));
        ClusterOwners(r(c),:) = [0,0];
        for k=2:length(cluster)
            [r,c] = find(ClusterOwners==common);
            order = [order cluster(r)];
            common = setdiff(ClusterOwners(r,:),ClusterOwners(r,c));
            ClusterOwners(r,:) = [0,0];
        end
        catch
        ClusterOwners = owners(Bloc(cluster,1),:);
        [blah,r] = max(ClusterMaxPos);
        [blah2,c] = max(blah);
        order = cluster(r(c));
        common = setdiff(ClusterOwners(r(c),:),ClusterOwners(r(c),c));
        ClusterOwners(r(c),:) = [0,0];
        for k=2:length(cluster)
            [r,c] = find(ClusterOwners==common);
            order = [order cluster(r)];
            common = setdiff(ClusterOwners(r,:),ClusterOwners(r,c));
            ClusterOwners(r,:) = [0,0];
        end
        order = fliplr(order);
        end
        Bexit(cluster) = Bexit(order);
        Bloc(cluster,:) = Bloc(order,:);
        Bexit(cluster) = round(Bexit(cluster)*10^3)/10^3 + (0:0.001:0.001*(length(cluster)-1));
    end
end

clear indy suspects cluster ClusterOwners ClusterMaxPos blah blah2 r c order common 
clear tempy own1 own2 only1 only2 yc y1 y2 xc x1 x2

Breaks = [Lf*ones(size(Rexit)) Texit zeros(size(Lexit)) Bexit; Rexit Wf*ones(size(Texit)) Lexit zeros(size(Bexit))];
Locs = [Rloc; Tloc; Lloc; Bloc];

if isempty(Breaks)
    ChannelWalls = cell(1,1);
    owners2 = 1;
    yr = (0:e:Wf); xr = Lf*ones(size(yw));
    xt = (Lf:-e:0); yt = Wf*ones(size(xt));
    yl = (Wf:-e:0); xl = zeros(size(yl));
    xb = (0:e:Lf); yb = zeros(size(xb));
    traj = [xr xt xl xb; yr yt yl yb];
    ChannelWalls{1,1} = traj;
else
    ChannelWalls = cell(length(Rexit)+length(Texit)+length(Lexit)+length(Bexit),1);
    owners2 = zeros(size(ChannelWalls));
    % create edge trajectories
    Breaks(:,end+1) = Breaks(:,1);
    Locs(end+1,:) = Locs(1,:);
    for n=1:size(Breaks,2)-1
        % x-coord the same
        if Breaks(1,n)==Breaks(1,n+1) && Breaks(1,n)==Lf && Breaks(2,n)<Breaks(2,n+1) % right edge, go up
            yw = (Breaks(2,n):e:Breaks(2,n+1));
            xw = Breaks(1,n)*ones(size(yw));
            traj = [xw; yw];
            ChannelWalls{n,1} = traj;
        elseif Breaks(1,n)==Breaks(1,n+1) && Breaks(1,n)==0 && Breaks(2,n)>Breaks(2,n+1) % left edge, go down
            yw = (Breaks(2,n):-e:Breaks(2,n+1));
            xw = Breaks(1,n)*ones(size(yw));
            traj = [xw; yw];
            ChannelWalls{n,1} = traj;           
        % y-coord the same
        elseif Breaks(2,n)==Breaks(2,n+1) && Breaks(2,n)==0 && Breaks(1,n)<Breaks(1,n+1) % bottom edge, go right
            xw = (Breaks(1,n):e:Breaks(1,n+1));
            yw = Breaks(2,n)*ones(size(xw));
            traj = [xw; yw];
            ChannelWalls{n,1} = traj;
        elseif Breaks(2,n)==Breaks(2,n+1) && Breaks(2,n)==Wf && Breaks(1,n)>Breaks(1,n+1) % top edge, go left
            xw = (Breaks(1,n):-e:Breaks(1,n+1));
            yw = Breaks(2,n)*ones(size(xw));
            traj = [xw; yw];
            ChannelWalls{n,1} = traj;
        else % turning corner(s)
            if Breaks(1,n)==Lf % from right edge
                if Breaks(1,n+1)==Lf % to right edge
                    yr = (Breaks(2,n):e:Wf); xr = Lf*ones(size(yr));
                    xt = (Lf:-e:0); yt = Wf*ones(size(xt));
                    yl = (Wf:-e:0); xl = zeros(size(yl));
                    xb = (0:e:Lf); yb = zeros(size(xb));
                    yr2 = (0:e:Breaks(2,n+1)); xr2 = Lf*ones(size(yr2));
                    traj = [xr xt xl xb xr2; yr yt yl yb yr2];
                    ChannelWalls{n,1} = traj;
                elseif Breaks(2,n+1)==Wf % to top edge
                    yr = (Breaks(2,n):e:Wf); xr = Lf*ones(size(yr));
                    xt = (Lf:-e:Breaks(1,n+1)); yt = Wf*ones(size(xt));
                    traj = [xr xt; yr yt];
                    ChannelWalls{n,1} = traj;
                elseif Breaks(1,n+1)==0 % to left edge
                    yr = (Breaks(2,n):e:Wf); xr = Lf*ones(size(yr));
                    xt = (Lf:-e:0); yt = Wf*ones(size(xt));
                    yl = (Wf:-e:Breaks(2,n+1)); xl = zeros(size(yl));
                    traj = [xr xt xl; yr yt yl];
                    ChannelWalls{n,1} = traj;
                elseif Breaks(2,n+1)==0 % to bottom edge
                    yr = (Breaks(2,n):e:Wf); xr = Lf*ones(size(yr));
                    xt = (Lf:-e:0); yt = Wf*ones(size(xt));
                    yl = (Wf:-e:0); xl = zeros(size(yl));
                    xb = (0:e:Breaks(1,n+1)); yb = zeros(size(xb));
                    traj = [xr xt xl xb; yr yt yl yb];
                    ChannelWalls{n,1} = traj;
                end
            elseif Breaks(2,n)==Wf % from top edge
                if Breaks(1,n+1)==Lf % to right edge
                    xt = (Breaks(1,n):-e:0); yt = Wf*ones(size(xt));
                    yl = (Wf:-e:0); xl = zeros(size(yl));
                    xb = (0:e:Lf); yb = zeros(size(xb));
                    yr = (0:e:Breaks(2,n+1)); xr = Lf*ones(size(yr));
                    traj = [xt xl xb xr; yt yl yb yr];
                    ChannelWalls{n,1} = traj;
                elseif Breaks(2,n+1)==Wf % to top edge
                    xt = (Breaks(1,n):-e:0); yt = Wf*ones(size(xt));
                    yl = (Wf:-e:0); xl = zeros(size(yl));
                    xb = (0:e:Lf); yb = zeros(size(xb));
                    yr = (0:e:Wf); xr = Lf*ones(size(yr));
                    xt2 = (Lf:-e:Breaks(1,n+1)); yt2 = Wf*ones(size(xt2));
                    traj = [xt xl xb xr xt2; yt yl yb yr yt2];
                    ChannelWalls{n,1} = traj;
                elseif Breaks(1,n+1)==0 % to left edge
                    xt = (Breaks(1,n):-e:0); yt = Wf*ones(size(xt));
                    yl = (Wf:-e:Breaks(2,n+1)); xl = zeros(size(yl));
                    traj = [xt xl; yt yl];
                    ChannelWalls{n,1} = traj;
                elseif Breaks(2,n+1)==0 % to bottom edge
                    xt = (Breaks(1,n):-e:0); yt = Wf*ones(size(xt));
                    yl = (Wf:-e:0); xl = zeros(size(yl));
                    xb = (0:e:Breaks(1,n+1)); yb = zeros(size(xb));
                    traj = [xt xl xb; yt yl yb];
                    ChannelWalls{n,1} = traj;
                end                
            elseif Breaks(1,n)==0 % from left edge
                if Breaks(1,n+1)==Lf % to right edge
                    yl = (Breaks(2,n):-e:0); xl = zeros(size(yl));
                    xb = (0:e:Lf); yb = zeros(size(xb));
                    yr = (0:e:Breaks(2,n+1)); xr = Lf*ones(size(yr));
                    traj = [xl xb xr; yl yb yr];
                    ChannelWalls{n,1} = traj;
                elseif Breaks(2,n+1)==Wf % to top edge
                    yl = (Breaks(2,n):-e:0); xl = zeros(size(yl));
                    xb = (0:e:Lf); yb = zeros(size(xb));
                    yr = (0:e:Wf); xr = Lf*ones(size(yr));
                    xt = (Lf:-e:Breaks(1,n+1)); yt = Wf*ones(size(xt));
                    traj = [xl xb xr xt; yl yb yr yt];
                    ChannelWalls{n,1} = traj;
                elseif Breaks(1,n+1)==0 % to left edge
                    yl = (Breaks(2,n):-e:0); xl = zeros(size(yl));
                    xb = (0:e:Lf); yb = zeros(size(xb));
                    yr = (0:e:Wf); xr = Lf*ones(size(yr));
                    xt = (Lf:-e:0); yt = Wf*ones(size(xt));
                    yl2 = (0:-e:Breaks(2,n+1)); xl2 = zeros(size(yl2));
                    traj = [xl xb xr xt xl2; yl yb yr yt yl2];
                    ChannelWalls{n,1} = traj;
                elseif Breaks(2,n+1)==0 % to bottom edge
                    yl = (Breaks(2,n):-e:0); xl = zeros(size(yl));
                    xb = (0:e:Breaks(1,n+1)); yb = zeros(size(xb));
                    traj = [xl xb; yl yb];
                    ChannelWalls{n,1} = traj;
                end
            elseif Breaks(2,n)==0 % from bottom edge
                if Breaks(1,n+1)==Lf % to right edge
                    xb = (Breaks(1,n):e:Lf); yb = zeros(size(xb));
                    yr = (0:e:Breaks(2,n+1)); xr = Lf*ones(size(yr));
                    traj = [xb xr; yb yr];
                    ChannelWalls{n,1} = traj;
                elseif Breaks(2,n+1)==Wf % to top edge
                    xb = (Breaks(1,n):e:Lf); yb = zeros(size(xb));
                    yr = (0:e:Wf); xr = Lf*ones(size(yr));
                    xt = (Lf:-e:Breaks(1,n+1)); yt = Wf*ones(size(xt));
                    traj = [xb xr xt; yb yr yt];
                    ChannelWalls{n,1} = traj;
                elseif Breaks(1,n+1)==0 % to left edge
                    xb = (Breaks(1,n):e:Lf); yb = zeros(size(xb));
                    yr = (0:e:Wf); xr = Lf*ones(size(yr));
                    xt = (Lf:-e:0); yt = Wf*ones(size(xt));
                    yl = (Wf:-e:Breaks(2,n+1)); xl = zeros(size(yl));
                    traj = [xb xr xt xl; yb yr yt yl];
                    ChannelWalls{n,1} = traj;
                elseif Breaks(2,n+1)==0 % to bottom edge
                    xb = (Breaks(1,n):e:Lf); yb = zeros(size(xb));
                    yr = (0:e:Wf); xr = Lf*ones(size(yr));
                    xt = (Lf:-e:0); yt = Wf*ones(size(xt));
                    yl = (Wf:-e:0); xl = zeros(size(yl));
                    xb2 = (0:e:Breaks(1,n+1)); yb2 = zeros(size(xb2));
                    traj = [xb xr xt xl xb2; yb yr yt yl yb2];
                    ChannelWalls{n,1} = traj;
                end
            end
        end
    end
    % identify owner max of edge domains
    loner = [];
    for n=1:size(Breaks,2)-1
        loc1 = Locs(n,:);
        loc2 = Locs(n+1,:);
        own1 = owners(loc1(1),:);
        own2 = owners(loc2(1),:);
        common = intersect(own1,own2);
        if length(common)==1
            owners2(n) = common;
        elseif loc1(1)==loc2(1) && loc1(2)~=loc2(2) && own1(1)==own2(1) && own1(2)==own2(2)
            loner = [loner; n own1];
        else
            if Breaks(1,n)==Lf
                init = 'R';
            elseif Breaks(1,n)==0
                init = 'L';
            elseif Breaks(2,n)==Wf
                init = 'T';
            elseif Breaks(2,n)==0
                init = 'B';
            end
            if Breaks(1,n+1)==Lf
                fin = 'R';
            elseif Breaks(1,n+1)==0
                fin = 'L';
            elseif Breaks(2,n+1)==Wf
                fin = 'T';
            elseif Breaks(2,n+1)==0
                fin = 'B';
            end
            trajI1 = C2{loc1(1),1};
            trajI2 = C2{loc1(1),2};
            trajF1 = C2{loc2(1),1};
            trajF2 = C2{loc2(1),2};
            
            if strcmp(init,'R')
                if trajI1(2,10)>saddles(loc1(1),2) && trajI2(2,10)<saddles(loc1(1),2)
                    passI = own1(1);
                elseif trajI2(2,10)>saddles(loc1(1),2) && trajI1(2,10)<saddles(loc1(1),2)
                    passI = own1(2);
                else
                    figure; text(1,1,'SOS!','FontSize',20);
                    error('test inconsistent')
                end
            elseif strcmp(init,'L')
                if trajI1(2,10)<saddles(loc1(1),2) && trajI2(2,10)>saddles(loc1(1),2)
                    passI = own1(1);
                elseif trajI2(2,10)<saddles(loc1(1),2) && trajI1(2,10)>saddles(loc1(1),2)
                    passI = own1(2);
                else
                    figure; text(1,1,'SOS!','FontSize',20);
                    error('test inconsistent')
                end
            elseif strcmp(init,'T')
                if trajI1(1,10)<saddles(loc1(1),1) && trajI2(1,10)>saddles(loc1(1),1)
                    passI = own1(1);
                elseif trajI2(1,10)<saddles(loc1(1),1) && trajI1(1,10)>saddles(loc1(1),1)
                    passI = own1(2);
                else
                    figure; text(1,1,'SOS!','FontSize',20);
                    error('test inconsistent')
                end
            elseif strcmp(init,'B')
                if trajI1(1,10)>saddles(loc1(1),1) && trajI2(1,10)<saddles(loc1(1),1)
                    passI = own1(1);
                elseif trajI2(1,10)>saddles(loc1(1),1) && trajI1(1,10)<saddles(loc1(1),1)
                    passI = own1(2);
                else
                    figure; text(1,1,'SOS!','FontSize',20);
                    error('test inconsistent')
                end
            end
            
            if strcmp(fin,'R')
                if trajF1(2,10)<saddles(loc2(1),2) && trajF2(2,10)>saddles(loc2(1),2)
                    passF = own2(1);
                elseif trajF2(2,10)<saddles(loc2(1),2) && trajF1(2,10)>saddles(loc2(1),2)
                    passF = own2(2);
                else
                    figure; text(1,1,'SOS!','FontSize',20);
                    error('test inconsistent')
                end
            elseif strcmp(fin,'L')
                if trajF1(2,10)>saddles(loc2(1),2) && trajF2(2,10)<saddles(loc2(1),2)
                    passF = own2(1);
                elseif trajF2(2,10)>saddles(loc2(1),2) && trajF1(2,10)<saddles(loc2(1),2)
                    passF = own2(2);
                else
                    figure; text(1,1,'SOS!','FontSize',20);
                    error('test inconsistent')
                end
            elseif strcmp(fin,'T')
                if trajF1(1,10)>saddles(loc2(1),1) && trajF2(1,10)<saddles(loc2(1),1)
                    passF = own2(1);
                elseif trajF2(1,10)>saddles(loc2(1),1) && trajF1(1,10)<saddles(loc2(1),1)
                    passF = own2(2);
                else
                    figure; text(1,1,'SOS!','FontSize',20);
                    error('test inconsistent')
                end
            elseif strcmp(fin,'B')
                if trajF1(1,10)<saddles(loc2(1),1) && trajF2(1,10)>saddles(loc2(1),1)
                    passF = own2(1);
                elseif trajF2(1,10)<saddles(loc2(1),1) && trajF1(1,10)>saddles(loc2(1),1)
                    passF = own2(2);
                else
                    figure; text(1,1,'SOS!','FontSize',20);
                    error('test inconsistent')
                end
            end
            
            if passI == passF
                owners2(n) = passI;
            else
                figure; text(1,1,'SOS!','FontSize',20);
                error('test inconsistent')
            end
            
        end
    end
    if ~isempty(loner)
        for n=1:size(loner,1)
            r1 = find(owners2==loner(n,2));
            r2 = find(owners2==loner(n,3));
            if isempty(r1) && ~isempty(r2)
                owners2(loner(n,1)) = loner(n,2);
            elseif isempty(r2) && ~isempty(r1)
                owners2(loner(n,1)) = loner(n,3);
            else
                figure; text(1,1,'SOS!','FontSize',20);
                error('*Now* Im confused!')
            end
        end
    end
end

clear Lexit Rexit Texit Bexit Lloc Rloc Tloc Bloc traj indy Breaks Locs
clear yr xt yl xb yw xw xr yt xl yb yr2 xr2 xt2 yt2 yl2 xl2 xb2 yb2 passF
clear loc1 loc2 own1 own2 common init fin trajI1 trajI2 trajF1 trajF2 passI
clear loner r1 r2

% step (4): integrate 1, u, u^2 over area of each domain. 
% ID Left/Right doms & calc openings.
% find max u on each domain.

area = zeros(length(maxima),1);
u_int = zeros(length(maxima),1);
u2_int = zeros(length(maxima),1);
Rdoms = [];
Ldoms = [];
Ropen = [];
Lopen = [];
MaxU = zeros(length(maxima),1);

for n=1:length(maxima)
    [r,c] = find(owners==n);
    walls = [];
    for k=1:length(r)
        traj = C{r(k),1};
        walls = [walls, traj];
        traj = C{r(k),2};
        walls = [walls, traj];
    end
    walls(3,:)=[];
    r = find(owners2==n);
    for k=1:length(r)
        traj = ChannelWalls{r(k)};
        walls = [walls, traj];
    end
    
    r = find(walls(1,:)==0);
    if ~isempty(r)
        Ldoms = [Ldoms n];
        Lopen = [Lopen e*(length(r)-1)];
    end
    r = find(walls(1,:)==Lf);
    if ~isempty(r)
        Rdoms = [Rdoms n];
        Ropen = [Ropen e*(length(r)-1)];
    end
    
    box = [round(min(walls(1,:))/e)*e round(max(walls(1,:))/e)*e round(min(walls(2,:))/e)*e round(max(walls(2,:))/e)*e];
    x = (box(1):e:box(2));
    y = (box(3):e:box(4));
    Umat = zeros(length(y),length(x));
    
    walls = round(walls/e)*e;
    IntBnds = zeros(length(y),2);
    for k=1:length(y)
        r = find(round(walls(2,:)/e)==round(y(k)/e));
        if ~isempty(r)
            blah = min(walls(1,r));
            IntBnds(k,1) = find(round(x/e)==round(blah/e));
            blah = max(walls(1,r));
            IntBnds(k,2) = find(round(x/e)==round(blah/e));
        end
    end
    
    [r,c] = find(IntBnds==0);
    for k=1:length(r)
        IntBnds(r(k),c(k)) = (IntBnds(r(k)-1,c(k)));
    end
    
    for k=1:size(keep,1)
        xk = (keep(k,1):e:keep(k,2));
        yk = (keep(k,3):e:keep(k,4));
        xc = intersect(round(xk/e),round(x/e))*e;
        yc = intersect(round(yk/e),round(y/e))*e;
        if ~isempty(xc) && ~isempty(yc)
            u = info{k,2};
            [Xc,Yc] = meshgrid(xc,yc);
            ue = feval(u,Xc,Yc);
            Umat(round((yc-box(3))/e+1),round((xc-box(1))/e+1)) = ue;
        end
    end
    
    I = zeros(size(Umat));
    Umat2 = zeros(size(Umat));
    for k=1:length(y)
        Umat2(k,IntBnds(k,1):IntBnds(k,2)) = Umat(k,IntBnds(k,1):IntBnds(k,2));
        I(k,IntBnds(k,1):IntBnds(k,2)) = 1;
    end
    
    MaxU(n) = max(max(Umat2));
    u_int(n) = trapz(x,trapz(y,Umat2,1));
    u2_int(n) = trapz(x,trapz(y,Umat2.^2,1));
    area(n) = trapz(x,trapz(y,I,1));
end

clear r c walls traj box x y Umat IntBnds blah xk yk xc yc u Xc Yc ue I Umat2

% step (5): calc on-site energies

localEs = u_int./u2_int;

% step (6): identify links btw nearest neighbours

owners = sort(owners,2,'ascend');
[nns,links,~] = unique(owners,'rows');

for k=1:size(nns,1)
    here=2;
    for n=links(k,1)+1:size(owners,1)
        if owners(n,:)==nns(k,:)
            links(k,here) = n;
            here = here + 1;
        end
    end
end

clear here

% step (7): calculate energy-dependant \rho_ij for nn's.

rhoij = zeros(size(nns,1),length(Evec));

for n=1:size(nns,1)
    paths = links(n,:);
    rm = find(paths==0);
    paths(rm) = [];
    for k=1:length(Evec)
        options = zeros(size(paths));
        for j=1:length(options)
            traj1 = C2{paths(j),1};
            traj2 = C2{paths(j),2};
            options(j) = e*trapz(real(sqrt(1./traj1(3,:)-Evec(k)))) + e*trapz(real(sqrt(1./traj2(3,:)-Evec(k))));
        end
        rhoij(n,k) = mean(options);
    end
end

clear paths rm traj1 traj2 options

% step (8): calculate norms for nn's.

Norms = zeros(size(maxima,1),length(Evec));
for n=1:size(maxima,1)
    [r,c] = find(nns==n);
    for j=1:length(r)
        c(j) = setdiff([1,2],c(j));
    end
    yada = zeros(1,length(Evec));
    for j=1:length(r)
        yada = yada + area(nns(r(j),c(j)))*exp(-2*real(rhoij(r(j),:)));
    end   
    Norms(n,:) = sqrt(u2_int(n) + yada*u_int(n)^2/area(n)^2);
end

% step (9): calculate nn tunneling rates.

tij = zeros(size(rhoij));
for n=1:size(nns,1)
    tij(n,:) = exp(-rhoij(n,:))./Norms(nns(n,1),:)./Norms(nns(n,2),:)*u_int(nns(n,1))*u_int(nns(n,2))*(localEs(nns(n,1))/area(nns(n,2)) + localEs(nns(n,2))/area(nns(n,1)));
end

% this bit saves the resulting workspace. I've commented it out for you as I've uploaded the processed workspaces as well, and this would overwrite it.
if 0
folder = strcat('L',int2str(L),'W',int2str(W),'ff',num2str(ff),'P');
folder = strrep(folder,'.','_');
cd SimCity; % choose folder
if ~exist(folder,'dir')
    mkdir(folder);
end
cd(folder);
save(strcat('wrkspc',int2str(count),'.mat'),'area','aTE','C','C2','ChannelWalls','count','e','E0','ell','Evec','extrema','ff','folder','hbar','info','kB','L','Ldoms','Lf',...
    'links','localEs','Lopen','m','maxima','MaxU','minima','nns','Norms','owners','owners2','posx','posy','Rdoms','rhoij','Ropen','saddles','T','tij','u2_int',...
    'u_int','W','Wf');
cd ../..
end
figure; text(1,1,'Done!','FontSize',20);
