% stability.m - generate various figures in this chapter
% RMM, 9 Oct 02
%
% 16 Apr 06, RMM: converted to standard plotting format
%
% 8 Jul 04, RMM: copied to AM04 directory
%   * converted print codes to generate EPS files
%   * reset some fonts and sizes to look right for the book
%   * Got rid of pauses
%
% 11 Oct 03, RMM: updated lecture examples
%   * updated code to make use of new version of phaseplot (with args)
%   * minimized the number of external functions used
%   * updated predator prey example
%
% Required files: oscillator.m, pendulum.m, saddle.m, predprey.m
%
% AM04 note: this code should all be pulled out into separate files so that
% it can be included on the web site as separate figures (I think).

aminit;					% initialize MATALB

%
% Stability definitions
%
% This set of plots illustrates the various types of equilibrium points.
%
fprintf(2, 'Stability:'); 

%
% Asy stable
%
fprintf(2, 'asystable...'); 
m = 1; b = 1; k = 1;			% default values
figure (1); clf; subplot(221);
amaxis([-1 1 -1 1]); set(gca, 'DataAspectRatio', [1 1 1]);
amphaseplot('oscillator', 'timepts', [0.3 1 2 3], nan, ...
  [-1,1; -0.3,1; 0,1; 0.25,1; 0.5,1; 0.7,1; 1,1; 1.3,1; ...
   1,-1; 0.3,-1; 0,-1; -0.25,-1; -0.5,-1; -0.7,-1; -1,-1; -1.3,-1], ...
  0:0.01:10, m, b, k);
hold on; plot([0], [0], 'k.', 'MarkerSize', AM_data_markersize*3);
% set(gca,'FontSize', 16); 
xlabel('{\itx}_1'); ylabel('{\itx}_2');

amprint('asystablepp.eps');

[t,x] = ode45('oscillator', [0 10], [1 0], [], m, b, k);
figure(2); clf; subplot(321);
plot(t, x(:,1), '-', t, x(:,2), '--', 'LineWidth', AM_data_linewidth);

% a=gca; set(a,'FontSize',20);
% set(gca,'DataAspectRatio',[2,1,1]);
amaxis;
xlabel('Time {\itt}'); ylabel('{\itx}_1, {\itx}_2');

legend('{\itx}_1', '{\itx}_2', 'Orientation', 'Horizontal'); 
lgh = legend; legend(lgh, 'boxoff');

amprint('asystablesim.eps');

%
% Unstable
%
% Note: parameter values need tweaking to generate plot from lecture
%
fprintf(2, 'saddle...'); 
figure(3); clf; subplot(221);
amaxis([-1 1 -1 1]); set(gca, 'DataAspectRatio', [1 1 1]);
amphaseplot('saddle', 'timepts', [0.2 0.5 0.8], nan, ...
  [-1, -1; 1, 1; ...
   -1, -0.95; -1, -0.9; -1, -0.8; -1, -0.6; -1, -0.4; -1, -0.2;...
   -0.95, -1; -0.9, -1; -0.8, -1; -0.6, -1; -0.4, -1; -0.2, -1;...
   1, 0.95; 1, 0.9; 1, 0.8; 1, 0.6; 1, 0.4; 1, 0.2;...
   0.95, 1; 0.9, 1; 0.8, 1; 0.6, 1; 0.4, 1; 0.2, 1;...
   -0.5, -0.45; -0.45, -0.5; 0.5, 0.45; 0.45, 0.5; ...
  -0.04, 0.04; 0.04, -0.04 ], 0:0.01:2);
hold on; plot([0], [0], 'k.', 'MarkerSize', AM_data_markersize*3);
% set(gca,'FontSize', 16); 
xlabel('{\itx}_1'); ylabel('{\itx}_2');

amprint('unstablepp.eps');

[t,x] = ode45('saddle', [0 10], [1 0]);
figure(4); clf; subplot(321);
plot(t, x(:,1), '-', t, x(:,2), '--', 'LineWidth', AM_data_linewidth); 
amaxis([0 3 -100 100]);

% a=gca; set(a,'FontSize',20);
% set(gca,'DataAspectRatio',[1,50,1]);
xlabel('Time {\itt}'); ylabel('{\itx}_1, {\itx}_2');
legend('{\itx}_1', '{\itx}_2'); 
lgh = legend; legend(lgh, 'boxoff');

amprint('unstablesim.eps');

%
% Stable isL
%
fprintf(2, 'isL...');
m = 1; b = 0; k = 1;			% zero damping
figure(5); clf; subplot(221);
amaxis([-1 1 -1 1]); set(gca, 'DataAspectRatio', [1 1 1]);
amphaseplot('oscillator', 'timepts', ...
  [pi/6 pi/3 pi/2 2*pi/3 5*pi/6 pi 7*pi/6 4*pi/3 9*pi/6 5*pi/3 11*pi/6 2*pi], nan, ...
  [0.2,0; 0.4,0; 0.6,0; 0.8,0; 1,0; 1.2,0; 1.4,0], ...
  0:0.01:20, m, b, k);
hold on; plot([0], [0], 'k.', 'MarkerSize', AM_data_markersize*3);
% set(gca,'FontSize', 16); 
xlabel('{\itx}_1'); ylabel('{\itx}_2');

amprint('stablepp.eps');

[t,x] = ode45('oscillator', [0 10], [1 0], [], m, b, k);
figure(6); clf; subplot(321); 
plot(t, x(:,1), '-', t, x(:,2), '--', 'LineWidth', AM_data_linewidth); 
% a=gca; set(a,'FontSize',20);
% set(gca,'DataAspectRatio',[2,1,1]);

amaxis([0 10 -2.5 2.5]);
xlabel('Time {\itt}'); ylabel('{\itx}_1, {\itx}_2');
legend('{\itx}_1', '{\itx}_2', 'Orientation', 'Horizontal'); 
lgh = legend; legend(lgh, 'boxoff');

amprint('stablesim.eps');
