% Plane Wave Decomposition Method for the stadium potential with sides of length
% L, where the radius of the circular ends is unity.
%
% This version attempts to normalize the wave function before calculating the tension.
% I was attempting to solve the problem that pwdm yields a tension as a function of k 
% that is quite rough, which makes it difficult to find the minima. There was no
% improvement in the "roughness", but I got better dynamic range in the tension, so
% may as well keep this feature.
%
% 2015-Nov-29  C. Kankelborg
%
function [tension, thetas, phases, Coeffs, x1, y1] = pwdm3(k, L, b)

% tension = The badness of this choice of k. Ideally, tension would approach zero for k equal to an eigenfrequency.
% k = wavenumber
% L = stadium side length (L=0 is the circular stadium).
% b = density of discretization (Recommend minimum 3; 12 is a reasonable number for non-circular stadia. See Li et al.)

% How many plane wave components should we have?
Circumference = 2*(L+pi);
lambda = 2*pi/k;
Nthetas = round(b*Circumference/lambda);

% Decide how many primary nodes there should be:
Mnodes = round(Nthetas*5/3); % Li et al. used M=N, but in a previous paper used M = (5/3)*N.

% Select Mnodes-1 boundary points ("primary nodes").
S0 = (1:(Mnodes-1))./(Mnodes-1); % Mnodes-1 arc positions along stadium boundary
[x0,y0] = boundarypoints(S0, L); % Boundary point Cartesian coordinates
% To this list of Mnodes-1 points, prepend one interior point that is unlikely to be a node.
% This is going to make renormalization necessary later (look for ##, below).
x0 = [1/pi, x0];
y0 = [1/e,  y0];

% Create vector of boundary values, i.e., desired wave function values at the Mnodes primary nodes
Bvalues = [1; zeros(Mnodes-1,1)]; % note the initial 1, which is at the sole interior point.

% Select 3*(Mnodes-1) additional boundary points ("secondary nodes").
S1 = [ (1:(Mnodes-1))-0.25, (1:(Mnodes-1))-0.5, (1:(Mnodes-1))-0.75 ]./(Mnodes-1); 
   % (Mnodes-1) Arc positions along stadium boundary, spaced uniformly between those in S0.
Mnodes1 = numel(S1);
[x1,y1] = boundarypoints(S1, L); % Boundary point Cartesian coordinates

% Select Nthetas values of propagation angle and phase
thetas = pi*(1:Nthetas)/Nthetas;  % Equally spaced planewave angles
phases = 2*pi*rand(1,Nthetas);   % Random phases

% Construct the matrix M connecting plane-wave amplitudes to boundary values.
for i0 = 1:Mnodes
   M(i0,1:Nthetas) = cos(k * (x0(i0)*cos(thetas) + y0(i0)*sin(thetas)) + phases);
endfor

% Invert the matrix equation M*Coeffs = Bvalues (might need SVD!). Thereby obtain coefficients for a 
% wavefunction that is zero at the primary nodes.
Coeffs = M \ Bvalues; % Note that if the matrix is singular, Octave computes a least squares solution.
% badness = sum( (M*Coeffs).^2 )/Nthetas %This should be small if all is working.

% Normalize the wave function to unit RMS. Until now, the amplitude has been controlled by letting
% a single arbitrary interior point have psi=1 (see above, ##).
% First, select random interior points:
Ninterior = 100;
[x2,y2] = interiorpoints(Ninterior, L);
% Evaluate the wave function at all x2,y2:
psi2 = wavefunction(x2, y2, k, thetas, phases, Coeffs);
% Now use x2,y2 to perform a MonteCarlo normalization integral:
normalization = sqrt(mean(psi2.*psi2));
Coeffs ./= normalization; % Apply normalization to the coefficients.

% Calculate the tension (sum of squares of the wave function evaluated at the secondary nodes).
tension=0;
for i1 = 1:Mnodes1
   tension += sum( Coeffs' .* cos(k * (x1(i1)*cos(thetas) + y1(i1)*sin(thetas)) + phases) ).^2;
endfor
tension /= Nthetas; % Not in Li's recipe, but normalizing seems right to me.

% end









