% Choose *approximately* N quasi-random interior points within the Buminovich stadium with 
% side length L, distributed with constant probability per unit area.
%
% 2015-Nov-29  C. Kankelborg
function [x,y] = interiorpoints(N,L)
% N = number of points to choose.
% L = Stadium side length (end cap radius = 1).

% Create a regular grid of points on a rectangle that circumscribes the stadium. 
A = 2*L+pi; % Area of the stadium.
delta = sqrt(A/N); % Grid spacing, sized to achieve the desired number of points.
x1d = -(1+L/2):delta:(1+L/2); % 1D array of x values
y1d = -1:delta:1;             % 1D array of y values
Nx = numel(x1d);
Ny = numel(y1d);
x2d = ones(Ny,1) * x1d;   % 2D array of x values
y2d = y1d' * ones(1,Nx);  % 2D array of y values

% Now move each x,y pair randomly within its (delta x delta) unit cell.
x2d += delta * rand(Ny,Nx); 
y2d += delta * rand(Ny,Nx); 

% Keep only those points inside the stadium. Due to the random displacements, the
% result is only approximately N.
%inside = find( (y2d<1) & (y2d>-1) & x2d < ( L/2 + sqrt(1 - y2d.^2) ) ); % & (x2d > ( -L/2 - sqrt(1 - y2d.^2) ))
   % The above expression should need the inequality that I have commented out. Why does it work without it?
% A more economical expression that works and actually makes sense:
inside = find( abs(y2d)<1 & abs(x2d) < ( L/2 + sqrt(1 - y2d.^2) ) );
x = x2d(inside);
y = y2d(inside);
