function p5b
% March 2016, March 2017, J. Gaspar

option = questdlg({'Consider friction/damping?', ...
    'hasFriction = nonzero beta', 'noFriction = zero beta'}, ...
    'Friction yes/no', ...
    'hasFriction', 'noFriction', 'hasFriction');
switch option
    case 'hasFriction',
        t3 % eigenvalues
        t4(1) % time solution (symbolic)
    case 'noFriction'
        t4(2) % time solution (symbolic) for
end


function A= mk_A( m, K, beta )
% make the 2x2 dynamics matrix
if nargin<1
    m=1; K=1.25; beta=3;
end
A= [0 1; -K/m -beta/m];


function t3
%lst= {[0 1; -1 0], [0 1; -1 -2], [0 1; -1 -4], mk_A};
lst= { mk_A }; % single example
Vsav= [];
Dsav= [];
for i=1:length(lst)
    A= lst{i};
    [V,D]= eig(A);

    % save eigenvalues
    Dsav= [Dsav [D(1,1); D(2,2)]];
    % enforce x coord of eigenvectors = 1 (NOT always possible!)
    V= V./repmat(V(1,:),2,1);
    Vinv= inv(V);
    Vsav= [Vsav V];
end
Dsav
Vsav


function t4(tstId)
if nargin<1, tstId= 1; end
if tstId==1
    % damping implies convergence to the origin [0; 0]
    m=1; K=1.25; beta=3; x0=[1; 1];
else
    % no damping implies oscilating forever
    m=1; K=1; beta=0; x0=[1; 1];
end
A= mk_A(m, K, beta)

syms s
x_of_t= simplify( ilaplace(inv(s*eye(2)-A))*x0 ); % varing constants formula
fprintf(1,'e^{At} = \n'); pretty(x_of_t)

% Show one trajectory
tRange= 0:.1:(5*2.5); % 5x the slowest time constant
figure(tstId); clf; hold on;
phase_plot(x_of_t, tRange)
xlabel('x_1'); ylabel('x_2'); drawnow

% show two axis
plot([0 0], [0 .5], 'k', 'linewidth',1)
plot([0 .5], [0 0], 'k', 'linewidth',1)

% Show many trajectories
[x10,x20]= meshgrid(-1:1, -1:1);
x0= [x10(:) x20(:)]';
sImAinv= ilaplace(inv(s*eye(2)-A));
for i=1:size(x0,2)
    x_of_t= sImAinv*x0(:,i);
    figure(tstId); phase_plot(x_of_t, tRange); drawnow
end

% draw the two eigenvectors
[V,D]= eig(A);
if max(abs(imag(D)))<eps
    % real eigenvalues imply easy to understand eigenvectors
    plot([0 V(1,1)], [0 V(2,1)], 'r', 'linewidth',2)
    plot([0 -V(1,1)], [0 -V(2,1)], 'r--', 'linewidth',2)
    plot([0 V(1,2)], [0 V(2,2)], 'g', 'linewidth',2)
    plot([0 -V(1,2)], [0 -V(2,2)], 'g--', 'linewidth',2)
end
xlabel('x_1'); ylabel('x_2')
grid on
axis equal


function phase_plot(x_of_t, tRange)
x= [];
for t= tRange
    xi= eval(x_of_t);
    x(:,end+1)= xi;
end
%plot(x') % time plot
plot(x(1,:),x(2,:),'b.-') % xy plot