 function [] = orbit_plot(z,DATA)
% orbit_plot solves the initial value
% problem using the supplied (z) values for
% the initial adjoints. It then generates
% some plots of the solution
%
% We solve an initial value problem forward in time
% for the combined 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) : range-angle
%
%   set initial conditions
%
      x0= [ 1.0 ;
            0.0 ;
            1.0 ;
           z(:) ;
            0  ];
   DATA.range = 1.;
  
   [T,X] = ode45('rhs',[0:DATA.tf/50:DATA.tf],x0,[],DATA);
%
%
  rp    = X(:,1);
  nup   = X(:,7);
  xp    = rp.*cos(nup);
  yp    = rp.*sin(nup);
%
% set the screen
%
  figure
  sa = max(rp) + .2;
  axis([-sa, sa, -sa, sa]);
%
% plot the probe's orbit
%
  plot(xp,yp,'r'); hold
%
%   Plot the Center Point
%   Add a label
%
  plot(0,0,'+r',0,0,'or')
  title('Solar-Sail: Min-time orbit')
%
% plot the final orbit
%
  ma = max(size(rp));
  ra = rp(ma);
  na = nup(ma);
  nu = [0: pi/50 : 2*pi];
  r = ones(size(nu));
  r = r.*ra;
  x = r .* cos(nu);
  y = r .* sin(nu);
  plot(x,y,'b')
%
% plot the initial orbit
%
  x =  cos(nu);
  y =  sin(nu);
  plot(x,y,'g')
%
% compute the thrust angles
%
 phi =  atan2(-X(:,5),-X(:,6));
 tmp = pi/2 + nup - phi;
 u = cos(tmp); v = sin(tmp);
 quiver(xp(1:3:ma),yp(1:3:ma),u(1:3:ma),v(1:3:ma))
 quiver(xp(1:3:ma),yp(1:3:ma),u(1:3:ma),v(1:3:ma),'o')
 axis('square'); 
 grid
%
% return
%
