  function P_loc = Orbit(t,A,P,phi,pcode)

% Pre: t is a row vector of "times" and (A,P,phi) defines the orbit
%
%     x(tau)      cos(phi)  sin(phi)     (A-P)/2 + ((A+P)/2)cos(tau)
%              =                      *
%     y(tau)     -sin(phi)  cos(phi)     sqrt(A*P)sin(tau)
%
%     pcode (optional)  is a valid plot string, e.g., '-', '.', '+', '*', etc.
%
% Post:
%      P_loc = [ x(t); y(t) ] and if nargin==5, this locus of points
%      is plotted.

   c = cos(t);
   s = sin(t);
   x = ((P-A)/2) + ((P+A)/2)*c;
   y = sqrt(A*P)*s;
   cphi = cos(phi);
   sphi = sin(phi);
   P_loc = [cphi sphi; -sphi cphi]*[x;y];
   if nargin == 5
      plot(P_loc(1,:),P_loc(2,:),pcode)
   end
