% Plane Wave Decomposition Method for the stadium potential with sides of length
% L, where the radius of the circular ends is unity.
%
% This version eschews random phases. I instead have phases of 0 and pi/2 at each 
% value of theta. This is an attempt to solve the problem that pwdm yields a
% tension as a function of k that is quite rough -- noise of order 30% for k~5.5.
% Unfortunately, my deterministic choice of phases does not improve things at all.
%
% CCK 2015-Nov-28 Based on pwdm.m. 

function [tension, thetas, phases, Coeffs, x1, y1] = pwdm2(k, L, b)

% tension = The badness of this choice of k. Ideally, tension would approach zero for k equal to an eigenfrequency.
% thetas = array of azimuthal angles used for the superposition of planewaves to make the wave function.
% phases = array of phases used for the superposition of planewaves to make the wave function.
% Coeffs = array of coefficients used for the superposition of planewaves to make the wave function.
% x1, y1 = array of secondary nodes around stadium boundary
% 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 = 2*round(b*Circumference/lambda/2);
   % I've uses Heller's recommended number of planewaves here, but I make sure
   % the number is even so that I can have half cosines (phase=0) and half sines (phase=pi/2).
% 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.
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/2)/(Nthetas/2);  % Equally spaced planewave angles.
thetas = [thetas,thetas];
   %Note the 2x redunancy.
phases = [zeros(1,Nthetas/2), ones(1,Nthetas/2)*pi/2];   
   % First half are cosines, second half are sines.

% 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.

% 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









