% This is a macro used to create the plots for numeric integration example
% in the discrete-time chapter

% Example 13.2


% Clear memory
clear all

% Define the input file names
exname = 'extdisc_ex2';    % the name of this macro
sysname = [exname, 'sys'];


% The simulation time (and step size for plotting)
Tm = 0.01;    % maximum sample rate
ninter = 10;  % number of intersample values to plot

t0 = 0;
tf = 5;
t = t0:(Tm/ninter):tf;

% Define the system size
n = 2;  % number of plant states

% Clear the figures
figure(1); clf


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is the continuous-time simulation

% used to easily pass simulation data
SysData = struct('n', n);

% Define the system initial conditions
x0 = [0.5; 0];

% Run the simulation
opts = '';
fprintf('Running case-study #%d\n', j);
[t, x] = ode23(sysname, t, x0, opts, SysData, 'deriv');

% Define the system outputs
y = zeros(length(t), 2);  % needs to be the size of the output vector
for i=1:length(t);
  tmp = feval(sysname, t(i), x(i,:), '', SysData, 'output');
  tmp = tmp(:);  % make sure it fits
  y(i, :) = tmp';  
end

% Plot the results 
h1a = plot(t, y(:,1)); hold on;
h2a = plot(t, y(:,2), '--'); hold on;

% Calculate the errors for a few different sample rates
Ts = Tm*[1, 5, 10];
ns = length(Ts);
figure(2); clf;
h = [];
st = [': '; '--'; '- '];
for j =1:ns

  % Calcualte the discrete-time sequence
  T = Ts(j);
  nits = tf/T + 1;
  xi = zeros(nits, n);
  w = zeros(nits, n);
  xi(1,:) = x0';

  % copy the samples
  ind_k = 1:nits; ind_k = ind_k(:);
  ind_t = (ind_k - 1)*ninter*(T/Tm) + 1;
  xi(ind_k, :) = y(ind_t, :);

  % determine the inter-sample error
  for k=1:(nits-1)
    x1 = xi(k, 1); x2 = xi(k, 2);
    d = 2*x1 - 0.1*x1^3;
    xi1 = x1 + T*(2*x2 + T*d)/2;
    xi2 = x2 + T*d;
    w(k, 1) = xi(k+1, 1) - xi1;
    w(k, 2) = xi(k+1, 2) - xi2;
  end

  figure(2);
  h = [stairs(T*(ind_k - 1), w(:,2), st(j, :)); h]; hold on;
end

figure(1);
h1b = stairs(T*(ind_k - 1), xi(:,1));
h2b = stairs(T*(ind_k - 1), xi(:,2), '--');
set([h1a; h1b; h2a; h2b], 'LineWidth', 1.5);
xlabel('t');

% print out the results
eval(sprintf('print -deps %s_a\n', exname));


figure(2);
axis([0 5 -0.15 0.2]);
set(h, 'LineWidth', 1.5);
xlabel('t');
ylabel('w_2', 'Rotation', 0);
legend(h, sprintf('T = %3.2f', Ts(3)), sprintf('T = %3.2f', Ts(2)),  sprintf('T = %3.2f', Ts(1)));

% print out the results
eval(sprintf('print -deps %s_b\n', exname));