% This is my attempt to implement the Plane Wave Decomposition Method of Heller,
% based on the description provided by Li, Robnik and Hu, Phys. Rev. E 57:4, 
% pp. 4096-4105 (1998). This is for the stadium potential with sides of length
% L, where the radius of the circular ends is unity.
%
% CCK 2015-Jun-16 Started coding.
% CCK 2015-Nov-05 Decided this would be a subroutine for calculating the tension. Will need
%     a main program to wrap this and search k-space for eigenfrequencies.
% CCK 2015-Dec-14 Added optional outputs to facilitate plotting the wave function.

function [tension, thetas, phases, Coeffs, x1, y1] = pwdm(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.
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.

% 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









