% Return a boundary point (x,y) corresponding to normalized arclength S on a
% stadium with unit radius and side length L.
function [x,y] = boundarypoints(S,L)
% [x,y] = cartesian coordinates of the boundary point(s) specified by S. If S is an
%     array, then x and y will each be row vectors with the same number of elements as S.
% S = normalized position along perimeter, CCW from the right end of the top straight 
%     side of the stadium (0 <= S <= 1). This may be an array.
% L = length of sides. L>0. Note that the radius of the end caps is 1.
%
% Example:
% octave:1> [x,y] = boundarypoints(0:0.01:1, 1.3);
% octave:2> plot(x,y)
% octave:3> axis equal
%
% CCK 2015-Nov-05

Circumference = 2*(L + pi);
S2 = S*Circumference; % I will work in non-normalized (but still dimensionless) coordinates.
for k=1:numel(S)
   if S2(k) < 0
      S(k)
      error("S<0 is not acceptable.")
   elseif S2(k) < L
      x(k) = L/2 - S2(k);
      y(k) = 1;
   elseif S2(k) < L+pi
      phi = S2(k) - L; % angular coordinate from top of left end cap
      x(k) = -L/2 -sin(phi);
      y(k) = cos(phi);
   elseif S2(k) < 2*L + pi
      x(k) = S2(k) - (1.5*L + pi);
      y(k) = -1;
   elseif S2(k) <= Circumference
      phi = S2(k) - 2*L - pi; % angular coordinate from bottom of right end cap
      x(k) = L/2 + sin(phi);
      y(k) = -cos(phi);
   else
      S(k)
      error("S>1 is not acceptable.")
   endif
end

