% plot_directionfield
%    helper function to plot the direction field for a system with a 
%    given coefficient matrix
% 
% 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_directionfield( varargin )

% required input are the axes and matrix
    figaxes = varargin{1};
    mat     = varargin{2};
% optional are the number of arrows in either grid direction and color
    if length( varargin ) > 2
        num_vec = varargin{3};
    else
        num_vec = 11;
    end
    if length( varargin ) > 3
        df_color = varargin{4};
    else
        df_color = 'k';
    end
% any further arguments are properties to pass to quiver
    if length( varargin ) > 4
        quiver_prop = varargin{5:length(varargin)};
    else
        quiver_prop = {'ShowArrowHead', 'on', 'AutoScaleFactor', .4};
    end
    
% get the axis limits
    axes( figaxes );
    axis_limits = axis();
    
% set up x and y vectors for the direction field
    x_vec = linspace( axis_limits(1), axis_limits(2), num_vec );
    y_vec = linspace( axis_limits(3), axis_limits(4), num_vec );

% define an x  and a y grid and the corresponding derivatives
    x_grid  = [ ];
    y_grid  = [ ];
    xp_grid = [ ];
    yp_grid = [ ];
    
    for j1=1:num_vec
        for j2=1:num_vec
            x_grid = [ x_grid x_vec(j1) ];
            y_grid = [ y_grid y_vec(j2) ];
            
            der = mat*[x_vec(j1); y_vec(j2)];
            xp_grid = [ xp_grid der(1) ];
            yp_grid = [ yp_grid der(2) ];
        end
    end
    
    dx_len = sqrt( xp_grid.^2 + yp_grid.^2 );
    
    quiver( figaxes, ...
            x_grid, y_grid, xp_grid./dx_len, yp_grid./dx_len, ...
            'Color', df_color, ...
            quiver_prop{1:length(quiver_prop)} );
end
