      function xdot =  rhs(t,x,FLAG,DATA)
%
%      compute the right-hand side of the state-adjoint system
%
%      x(1) : [r] radial - position
%      x(2) : [u] radial - velocity
%      x(3) : [v] transverse - velocity
%      x(4) : lambda-r
%      x(5) : lambda-u
%      x(6) : lambda-v
%      x(7) : angular position - only needed for plotting purposes
%
%
%  compute the thrust steering angle from the adjoints
%
      phi =  atan2(-x(5),-x(6));
      accm = DATA.acc / (1 - DATA.beta * t);

      xdot(1) = x(2);
      xdot(2) = (x(3)*x(3) - 1.0/x(1)) / x(1) + accm * sin(phi);
      xdot(3) = -x(2) * x(3) / x(1) + accm * cos(phi);
      xdot(4) = (x(5) * (x(3)*x(3) - 2.0/x(1)) - x(6)*x(2)*x(3))  / (x(1) *x(1));
      xdot(5) = x(6) * x(3) / x(1) - x(4);
      xdot(6) = (x(6) * x(2) - 2.0*x(5)*x(3)) / x(1);
%
  if DATA.range > 0
       xdot(7) = x(3) / x(1);
   end % if DATA.range
%
  xdot = xdot(:); % make it a column
%
% return
%
