% Try plotting an eigenmode using pwdm4.m, gridpoints.m, and wavefunction.m 


% First explore the region around some particular mode.
L=1
k = 18.92:0.002:18.96;
tension=0; %initialize to prevent confusion
start_time = datestr(now)
Nk = numel(k)
for ii = 1:Nk
   ii
   tension(ii) = pwdm4(k(ii),L,12);
endfor

figure(1)
hold off
plot(k,log10(tension),'r')


% Now find the lowest tension of the k-values sampled so far.
[tensionbest, ibest] = min(tension)
kbest = k(ibest)

% Calculate and plot the wave function.
dx = 0.01;
[tensionbest, thetas, phases, Coeffs, xbound, ybound] = pwdm4(kbest,L,12);
[x2d,y2d,x1d,y1d,mask] = gridpoints(dx,L);
figure(2)
hold off
psi = wavefunction(x2d, y2d, kbest, thetas, phases, Coeffs);
hold off
imagesc(x1d,y1d,psi.*mask)
hold on
axis equal
plot(xbound,ybound,'k.')
axis equal