% plot_solutiontraj
%    helper function to plot the solution trajectories in the phase
%    plane for a system with a given coefficient matrix
%    the input figaxes must have hold on set
% 
% by Gavin LaRose, 10/31/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_solutiontraj( varargin )

% required input are the axes and matrix
    figaxes = varargin{1};
    mat     = varargin{2};
% optional are the line style, initial conditions to use, and other
% properties to pass to the plot routine 
    if length( varargin ) > 2
        style = varargin{3};
    else
        style = '-b';
    end
    if length( varargin ) > 3
        iclist = varargin{4};
    else
        iclist = [ ];
    end
    if length( varargin ) > 4
        plot_prop = varargin{5:length(varargin)};
    else
        plot_prop = {'LineWidth', 2};
    end

    axes( figaxes );
    axis_limits = axis();
    x_step = (axis_limits(2) - axis_limits(1))/6;
    y_step = (axis_limits(4) - axis_limits(3))/6;
    
% find the eigenvalues and vectors
    [ evecs, lambdas ] = eig( mat );
    
    if length(iclist) == 0
        if real(lambdas(1,1)) < 0 | real(lambdas(2,2)) < 0
            xlist = axis_limits(1):x_step:axis_limits(2);
            ylist = axis_limits(3):y_step:axis_limits(4);

            for j1 = 1:length(xlist)
                iclist = vertcat( iclist, [ xlist(j1), axis_limits(4) ] );
                iclist = vertcat( iclist, [ xlist(j1), axis_limits(3) ] );
            end
            for j2 = 2:length(ylist)-1
                iclist = vertcat( iclist, [ axis_limits(1) ylist(j2) ] );
                iclist = vertcat( iclist, [ axis_limits(2) ylist(j2) ] );
            end
        elseif real(lambdas(1,1)) == 0 & real(lambdas(2,2)) == 0
            xlist = ((axis_limits(1)+axis_limits(2))/2):x_step: ...
                    axis_limits(2);
            ymid  = (axis_limits(3) + axis_limits(4))/2;
            for j1 = 1:length(xlist)
                iclist = vertcat( iclist, [ xlist(j1), ymid ] );
            end
        else
            for j1 = -1:1
                for j2 = -1:1
                    iclist = vertcat( iclist, [ .01*j1 .01*j2 ] );
                end
            end
        end                
    end

% generate and plot solutions
    for j1 = 1:length(iclist)
        [t_sol, x_sol] = ode45( @(t,x) [ mat*x ], [ 0 20 ], iclist(j1,:) );
        plot( figaxes, x_sol(:,1), x_sol(:,2), style, ...
              plot_prop{1:length(plot_prop)} );
    end
end

            % vec1 = evecs(1:2,1);
            % vec2 = evecs(1:2,2);
            % nvec = [ (arrow_length/norm(vec1))*vec1'; ...
            %          (arrow_length/norm(vec2))*vec2' ];
            % nvec
            % for j1=1:2
            %     if lambdas(j1,j1) > 0
            %         iclist = vertcat( iclist, ...
            %                           [ .01*nvec(j1,:); -.01*nvec(j1,:);
            %                             .01*nvec(j1,:) + [ 0 .01 ];
            %                             .01*nvec(j1,:) - [ 0 .01 ];
            %                             -.01*nvec(j1,:) + [ 0 .01 ];
            %                             -.01*nvec(j1,:) - [ 0 .01 ] ] );
            %     else
            %         iclist = vertcat( iclist, ...
            %                           [ axis_limits(1)*nvec(j1,:); ...
            %                             axis_limits(2)*nvec(j1,:); ...
            %                             .75*axis_limits(1)*nvec(j1,:); ...
            %                             1.25*axis_limits(2)*nvec(j1,:); ...
            %                             .75*axis_limits(2)*nvec(j1,:); ...
            %                             1.25*axis_limits(2)*nvec(j1,:) ] );
            %     end
            % end
