% Run directly after DiagH.m or DiagHwithWE.m. Do not clear workspace in between.
for n=1:many
evec = reshape(V(:,jones(n)),[Nx+2,Ny+2]).';
n % Print out mode number
E(n) % And its energy
figure; % Choose one of the following two plotting options:
%surface(x,y,log(abs(evec))); colormap('jet'); shading interp; grid on; hold on; % Plot absolute value on a log scale
surface(x,y,abs(evec)); colormap('jet'); shading interp; grid on; hold on; % Plot absolute value on a linear scale

if 0 % This allows you to plot the valley network on top. You will need to run ValleyNetwork.m FOR THE SAME NOISE REALISATION THAT YOU HAVE JUST DIAGONALISED *FIRST*, save the cell array 'C' and the matrix 'maxima' as a mat file, and load it into the workspace before running this bit.
for k=1:length(C)
    traj = C{k,1};
    c = find(traj(3,:)<1/E(n));
    plot3(traj(1,c),traj(2,c),traj(3,c),'k.');
    traj = C{k,2};
    c = find(traj(3,:)<1/E(n));
    plot3(traj(1,c),traj(2,c),traj(3,c),'w.');
end
% These few lines add domain labels. 
for k=1:size(maxima,1)
    plot(maxima(k,1),maxima(k,2),'ko');
    text(maxima(k,1),maxima(k,2),strcat('\leftarrow ',int2str(k)),'FontSize',20);
end
end
pause; close all; % Hit any key to go to the next pic.
end

