% rocketspring.m - Script plots some trajectories

% (c) Philip D Loewen, 2012-10-23


Tfinal = 3.75*pi;   % Users can change this. Have fun!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure('Name','Rocket-Spring Trajectories on [0,T]')

title('Rocket-Spring Trajectories on [0,T]');
hold on;
plot(128*[-1,1],[0,0],'k');     % Draw the x1-axis
plot([0,0],128*[-1,1],'k');     % Draw the x2-axis
xlabel('x_1');
ylabel('x_2');

options = odeset('RelTol',1.e-6,'AbsTol',1.e-8);

for theta=linspace(0,min([pi,Tfinal]),17)

  % Define the RHS for the controlled system dynamics, depending on theta
  xdot = @(t,xx) [xx(2);-xx(1)+sign(sin(t-theta))];

  % Solve the controlled dynamics
  [ttt,xxx] = ode45(xdot,[0,Tfinal],[0;0]);
  N = length(ttt);               % Note number of t-values returned

  % Draw the trajectories in solid red, add heavy dots on their final points.
  plot(xxx(:,1),xxx(:,2),'r');
  plot(xxx(N,1),xxx(N,2),'r.','MarkerSize',14);
  
  % Cheat: symmetry gives mirror-image trajectories. Make those ones green.
  plot(-xxx(:,1),-xxx(:,2),'g');
  plot(-xxx(N,1),-xxx(N,2),'g.','MarkerSize',14);
  
  end;

axis equal;
XmaxEst = Tfinal*2/pi;
set(gca,'XLim',1.1*XmaxEst*[-1,1],'YLim',1.1*XmaxEst*[-1,1]);

disp('Copy-paste these commands into the command window to produce printable output:');
disp(sprintf('figure(%d);',gcf));
disp('set(gcf,''PaperPosition'',[1,1,6,6])');
disp('set(gcf,''Renderer'',''painters'')');
disp('set(gcf,''RendererMode'',''manual'')');
disp('print -r600 -depsc rocketspring3pi2.eps');