% Survey k-space for modes of the Bunimovich stadium.
% 2015-Dec-16 CCK
% 2015-Dec-17 CCK fixed estimate of density of states.

function [tension, k] = survey(k1, k2, L, safety, outfile)
% k1 = initial k of survey
% k2 = final k of survey
% L = side length of stadium (end caps have unit radius)
% safety = safety factor for survey resolution. I naively estimate mean density
%          of modes in 2D k-space is As/(pi^2), where As = 2*L+pi is stadium 
%          geometrical area. Note that only the first quadrant of k-space counts.
%          The area that first quadrant of k-space, within a circle of radius k 
%          is pi*k^2/4. So the mean interval between modes is dk = 2*pi/((2*L+pi)*k). 
%          The safety factor, typically greater than unity, controls the stepping
%          of the survey in k: delta_k = 2*pi/((2*L+pi)*k*safety). Default: safety=40
% outfile = where to save the survey data. Default = 'survey.mat'
% tension = array of tension values for all tested values of k.
% k = array of all tested k values.

if !exist('safety','var') 
   safety = 40  % Default value
end

if !exist('outfile','var')
   outfile = 'survey.mat'  % Default output filename
end

threshold = 1e-3; % tension threshold below which roots of d(tension)/dk are considered significant.

n=1; % k counter
m=1; % mode counter
k(n) = k1;
tension(n) = pwdm4(k(n),L);

%Prepare for plotting:
dx = 0.01
[x2d,y2d,x1d,y1d,mask] = gridpoints(dx,L);

do % matlab would use repeat rather than do.
   n=n+1
   k(n) = k(n-1) + 2*pi/((2*L+pi)*k(n-1)*safety);
   tension(n) = pwdm4(k(n),L);
   if (n>2) && ( tension(n-1) < min([threshold, tension(n-2), tension(n)]) )
      % Print some diagnostics
      discovered_mode = m
      mode_frequency = k(n-1)
      
      % Save relevant parameters in mode structure
      mode(m).k = k(n-1);
      [mode(m).tension, mode(m).thetas, mode(m).phases, mode(m).Coeffs, xbound, ybound] = pwdm4(k(n-1),L);
      
      % Plot the wave function
      psi = wavefunction(x2d, y2d, mode(m).k, mode(m).thetas, mode(m).phases, mode(m).Coeffs);
      hold off
      imagesc(x1d,y1d,psi.*mask)
      hold on
      plot(xbound,ybound,'k.')
      axis equal      
      pause(0); % Argument is in seconds; if absent, then wait for a key to be pressed
      save(outfile, "L", "k", "tension", "mode") % Save progress up until now.
      
      m=m+1;
   end
until k(n) >= k2

