% NONLINEAR PHENOMENA, Spring 2013, UCSB (ECE 183 / ME 169 / PHYS 106)
% Katie Byl
% Lecture 8 (April 25, 2013)
%
% Code to demonstrate LINEAR 2D system trajectories in Phase Space.

% For best results, just play with the next two lines, only:
show_vector_parts = 0; % A flag to illustrate into eigen-motions
b=2.5; % "b=2" is critically damped. b<2 oscillates. b>2 has two, real roots

% (You should need to play with anything else below.)
m=1;
k=1;
A = [0 1; -k/m -b/m]  % State space matrix for "dx/dt = A*x"
[v,lambda] = eig(A)   % eigenvectors (v) and eigenvalues (diag of lambda)

% Plot the flow, in 2D Phase Space:
kv=.07; % scaling for "flow" vectors
figure(2); clf
[x,dx] = meshgrid(-10:1:10,-10:1:10);
plot([x(:) x(:)+dx(:)*kv]', [dx(:) dx(:)-b*dx(:)*kv-k*x(:)*kv]', 'r-',...
    'LineWidth',1)
hold on; axis equal; grid on
plot([x(:) ]', [dx(:) ]', 'k.','MarkerSize',6)
set(gca,'FontSize',20,'XTick',-20:5:20,'YTick',-20:5:20)

% Below, a complicated way to give a nice title:
strval = ['m = ',num2str(m), ', b = ', num2str(b), ', k = ', num2str(k),...
    '. \lambda_1 = ',num2str(lambda(1,1),'%.2f'),...
    ', \lambda_2 = ',num2str(lambda(2,2),'%.2f')];
title(strval(:)','FontSize',20)

% We will only plot the eigenvectors when they are REAL-valued:
if isreal(v)
    %input('Hit enter to plot EIGENVECTORS:')
    fprintf('Plotting EIGENVECTORS (in green)\n')
    for vi = 1:2
        plot(v(1,vi)*[-10 10],v(2,vi)*[-10 10],'g-','LineWidth',3,...
            'Color',[0 .7 0]);
    end
    text(-11,1.5,'eigenvectors','HorizontalAlignment','right',...
        'FontSize',20)
    plot([-15 -11],2.5+[0 0],'k-','Color',[0 .7 0],'LineWidth',2);
end
% plot NULLCLINES:
n1 = [1;0];
n2 = [-b;k];
n2 = (1/norm(n2))*n2;
bluegrey = [.4 .4 .7];
plot(n1(1)*[-10 10],n1(2)*[-10 10],'b--','LineWidth',2,...
    'Color',bluegrey);
plot(n2(1)*[-10 10],n2(2)*[-10 10],'r--','LineWidth',2,...
    'Color',bluegrey);
text(-11,4.5,'nullclines','HorizontalAlignment','right',...
    'FontSize',20)
plot([-15 -11],5.5+[0 0],'k--','Color',bluegrey,'LineWidth',2);


% Get ready to take mouse clicks, to plot trajectories in phase space:
t = 0:.05:30;
fprintf('Click to select at INITIAL CONDITION...\n');
fprintf('(Click anywhere where x>10 to STOP.)\n');
gg=ginput(1);

% Below, we are waiting for a mouse click with x value > 10 (to stop)
while gg(1)<10
    x0 = gg(:); % x0 = position AND velocity at t=0
    
    if (b/m)==2*sqrt(k/m) % CRITICALLY DAMPED (special case solution)
        
        % x(t) = (c1+c2*t)*e^(lambda*t), where both lambdas are the same
        c1 = x0(1);
        c2 = x0(2) - lambda(1,1)*c1;
        xt = (c1+c2*t).*exp(lambda(1,1)*t);
        vt = (c2+lambda(1,1)*(c1+c2*t)).*(exp(lambda(1,1)*t));
        x1 = c1*exp(lambda(1,1)*t);
        x2 = c2*t.*exp(lambda(1,1)*t);
        v1 = (c1*lambda(1,1))*exp(lambda(1,1)*t);
        v2 = (c2*(1+t)).*exp(lambda(1,1)*t);
        if show_vector_parts
            p1=plot(x1,v1,'b.-')
            p2=plot(x2,v2,'m.-')
        end
        
    else % Either underdamped (imag or complex) or overdamped (real roots)
        
        c = v\x0 % Find constants to match initial conditions
        
        % "real" used below to take care of numerical round-off errors:
        xt = real(c(1)*v(1,1)*exp(lambda(1,1)*t) + c(2)*v(1,2)*exp(lambda(2,2)*t));
        vt = real(c(1)*v(2,1)*exp(lambda(1,1)*t) + c(2)*v(2,2)*exp(lambda(2,2)*t));
        
        if (b/m)>2*sqrt(k/m) % for overdamped case, only
            
            % Break into eigenvector components (two decaying exponentials)
            x1 = real(c(1)*v(1,1)*exp(lambda(1,1)*t));
            x2 = real(c(2)*v(1,2)*exp(lambda(2,2)*t));
            v1 = real(c(1)*v(2,1)*exp(lambda(1,1)*t));
            v2 = real(c(2)*v(2,2)*exp(lambda(2,2)*t));
            
            if show_vector_parts % Show points on the green vectors
                p1=plot(x1,v1,'b.-')
                p2=plot(x2,v2,'m.-')
            end
        end
    end
    ui = [1:10]; % an index to TIME, for plotting
    p = plot(xt(ui),vt(ui),'k-','LineWidth',2);
    
    for n=20:50:length(t)
        ui = [1:n];
        % Below, we "reuse" the same plotting object, with new points:
        set(p,'XData',xt(ui),'YData',vt(ui));
        drawnow
    end
    gg = ginput(1); % Get a NEW initial condition. exits if x>10.
    
    if (b/m)>=2*sqrt(k/m) % for critically or overdamped cases, only
        
        if show_vector_parts % Set to a NEUTRAL color...
            set(p1,'Color',[.5 .5 .5]);
            set(p2,'Color',[.4 .4 .4]);
        end
    end
end

