% plot_localtraj
%    helper function to plot linearized solution trajectories in the
%    phase, near an equilibrium point, for a system with a given
%    coefficient matrix the input figaxes must have hold on set
% 
% by Gavin LaRose, 11/1/2012
% 
% (c)2012 Regents of the University of Michigan
% This work is licensed under the
% Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
% To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/3.0/
% or send a letter to
% Creative Commons,
% 444 Castro Street, Suite 900,
% Mountain View, CA, 94041, USA.
%
function plot_localtraj( varargin )

% required input are the axes, equilibrium point, and matrix
    figaxes = varargin{1};
    eq_pt   = varargin{2};
    mat     = varargin{3};

% optional are the radius around the point to which to graph solutions,
% the line style, initial conditions to use, and other properties to
% pass to the plot routine
    if length( varargin ) > 3
        traj_radius = varargin{4};
    else
        traj_radius = 0.25;
    end
    if length( varargin ) > 4
        style = varargin{5};
    else
        style = '-b';
    end
    if length( varargin ) > 5
        iclist = varargin{6};
    else
        iclist = [ ];
    end
    if length( varargin ) > 6
        plot_prop = varargin{7:length(varargin)};
    else
        plot_prop = {'LineWidth', 2};
    end

    axes( figaxes );
    
% find the eigenvalues and vectors
    [ evecs, lambdas ] = eig( mat );
    
    if length(iclist) == 0
        if real(lambdas(1,1)) < 0 | real(lambdas(2,2)) < 0
            % start with 12 points around a circle
            xlist = zeros(1, 12);
            ylist = zeros(1, 12);
            theta_step = pi/6;
            for j1 = 1:12
                iclist = vertcat( iclist, [ traj_radius*cos(j1*theta_step), ...
                                            traj_radius*sin(j1*theta_step) ] );
            end
        elseif real(lambdas(1,1)) == 0 & real(lambdas(2,2)) == 0
            for j1 = 1:6
                iclist = vertcat( iclist, [ j1*(traj_radius)/6, 0 ] );
            end
        else
            % start with points around and close to the point
            for j1 = -1:1
                for j2 = -1:1
                    iclist = vertcat( iclist, [ .01*j1 .01*j2 ] );
                end
            end
        end                
    end
    
    trsq = traj_radius^2;
% generate and plot solutions
    for j1 = 1:length(iclist)
        [t_sol, x_sol] = ode45( @(t,x) [ mat*x ], [ 0 20 ], iclist(j1,:) );
        u_plot = [ ];
        v_plot = [ ];
        for j2 = 1:length(x_sol(:,1))
            if x_sol(j2,1)^2 + x_sol(j2,2)^2 <= trsq
                u_plot = [ u_plot x_sol(j2,1)+eq_pt(1) ];
                v_plot = [ v_plot x_sol(j2,2)+eq_pt(2) ];
            else
                if length(u_plot) > 0
                    plot( figaxes, u_plot, v_plot, ...
                          style, plot_prop{1:length(plot_prop)} );
                    u_plot = [ ];
                    v_plot = [ ];
                end
            end
        end
        if length(u_plot) > 0
            plot( figaxes, u_plot, v_plot, ...
                  style, plot_prop{1:length(plot_prop)} );
        end
    end
end
