% Review data generated by survey.m
clear  % Get rid of any variables laying around.

load("survey/survey43.3-50.mat");
Nmodes = numel(mode)
kmin = min(k)
kmax = max(k)
Area = 2*L + pi

% Cumulative distribution of modes
Cum = 1:Nmodes;
figure(3)
hold off
plot([mode.k],log10(Cum),'b')
CumWeyl = (Area/(4*pi))*([mode.k].^2 - kmin.^2);
hold on
plot([mode.k],log10(CumWeyl),'r')

% Test symmetry of modes
Npoints = 400;
[xin,yin] = interiorpoints(Npoints,L);
for m=1:Nmodes
   psi   = wavefunction( xin,  yin, mode(m).k, mode(m).thetas, mode(m).phases, mode(m).Coeffs);
   psi_x = wavefunction(-xin,  yin, mode(m).k, mode(m).thetas, mode(m).phases, mode(m).Coeffs);
   psi_y = wavefunction( xin, -yin, mode(m).k, mode(m).thetas, mode(m).phases, mode(m).Coeffs);
   x_asymmetry = min( [mean( (psi-psi_x).^2 ), mean( (psi+psi_x).^2 )] );
   y_asymmetry = min( [mean( (psi-psi_y).^2 ), mean( (psi+psi_y).^2 )] );
   asymmetry(m) = max( [x_asymmetry, y_asymmetry] );
end
figure(4)
hold off
hist(asymmetry,100)
axis([0,max(asymmetry),-10,50])


% Inspect the "worst" mode!
[a_worst, m_worst] = max(asymmetry);

% Calculate and plot the wave function.
dx = 0.01;
[x2d,y2d,x1d,y1d,mask] = gridpoints(dx,L);
figure(5)
hold off
psi_worst = wavefunction(x2d, y2d, mode(m_worst).k, mode(m_worst).thetas, mode(m_worst).phases, mode(m_worst).Coeffs);
hold off
imagesc(x1d,y1d,psi_worst.*mask)
hold on
axis equal

% Resolve the peak.

figure(6)
hold off
iworst = find(k == mode(m_worst).k)
k1 = k(iworst-1)
k2 = k(iworst+1)
safety=160
outfile='/tmp/surveyreview.mat' % Not really looking to save this.
[tension2, k2] = survey(k1, k2, L, safety, outfile);

figure(7)
hold off
plot(k, log10(tension),'r')
axis([k(iworst-4),k(iworst+4)])
hold on
plot(k2, log10(tension2),'b')
