% a demonstration that considers the systems
%    x' = [ 1  9/4; 0  -2 ] x
% and
%    x' = [ -9/(4c - 9)  -(4c - 9)/4;  4c/(4c - 9)  0 ] x
% which are the nondimensionalized, linearized systems about the
% equilibrium solutions (0,0) and (c, 4c/(4c-9)) of the
% nonlinear model for TB infected australian opossums,
%    S' = (a - b) S + a I - beta I S
%    I' = -(alpha + b) I + beta I S,
% where S = the population of susceptible opossums and I = the
% population of infected opossums.  this model assumes that the
% disease is fatal (no members of the infected group reenter the
% population as recovered individuals).
%
% the phase plane is graphed, showing the linear behavior near each of
% the equlibrium solutions along with nonlinear trajectories.
%
% configuration variables are
%    show_point = 0 or 1 (0 = the zero equilibrium, 1 = the nonzero)
%    c_value    = the value of c (significant for the nonzero equilibrium)
% interesting values of c are c > 9/4 (2.25) (so that it is also has a
% positive y-coordinate), which give a stable spiral sink.  if c = 3, 
% we have 
%    lambda = -3/2 +/- i sqrt(3)/2, with
%    v = [ -3 +/- sqrt(3); 8 ]
%
% 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 Aussie_Opossum_Phase_Plane()
%-------------------------------------------------------------------------------
% configuration variables

% value of c
    c_value = 3;
    
% should we show the direction field?
    show_direction_field = 1;
% number of vectors to graph in each direction
    num_vec = 15;
    
% should we put a pause between graphing the direction field and linear
% trajectories, and graphing the nonlinear trajectories?
    include_pause = 1;
    
% the distance from the equilibrium solutions to plot the trajectories
% for the linearized systems
    linear_radius = 0.25;
    
% colors
    ltraj_color = 'b';  % the color for linear trajectories    
    ntraj_color = 'm';  % and for the nonlinear trajectories

    df_color    = 'k';  % the color of the direction field, if any
    
    background_grey = [ 0.6 0.6 0.6 ];

%-------------------------------------------------------------------------------
% initialization

% useful values
    r = 9/4;
    c = c_value;
    cfac = 4*c - 9;

% the non-zero equilibrium point
    nz_equil = [ c, 4*c/(4*c - 9) ];

% set up the figure window
    scrsz    = get(0,'ScreenSize');
    scrsz(3) = min(1500, scrsz(3));
    scrsz(4) = min(900, scrsz(4));
    
    fig1 = figure( 'OuterPosition', ...
                   [50 scrsz(4)/2 scrsz(3)/1.5 scrsz(3)/2.5], ...
                   'Name', 'Phase Plane' );
    figaxes = axes( 'Parent', fig1, 'Color', background_grey );
    hold(figaxes);
    
% put in some axes
    x_max = c + 1;
    y_max = round(nz_equil(2)) + 1;
    plot( figaxes, [-.5 x_max], [0 0], ':k', ...
          [0 0], [-.5 y_max], ':k', 'LineWidth', 2 );

    xlabel( figaxes, 'x', 'FontSize', 15 );
    ylabel( figaxes, 'y', 'FontSize', 15 );
    title( ['Phase Portrait for Australian Opossum System (c = ' ...
            num2str(c) ')'], 'FontSize', 15 );

    axis( [ -.5 x_max -.5 y_max ] );

%-------------------------------------------------------------------------------
% graph the trajectories for the linearized systems

    if show_direction_field
	% set up x and y vectors for the direction field
        x_vec = linspace( -0.5, x_max, num_vec );
        y_vec = linspace( -0.5, y_max, 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 = [ x_vec(j1) + r*y_vec(j2) - x_vec(j1)*y_vec(j2); ...
                        -c*y_vec(j2) + 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, 'ShowArrowHead', 'on', ...
                'AutoScaleFactor', .3 );
    end

% next plot the local trajectories near each equilibrium point
    plot_localtraj( figaxes, [0 0], [ 1 r; 0 -c ], linear_radius, ...
                    [ '-' ltraj_color ] );
    plot_localtraj( figaxes, nz_equil, [  -9/cfac -cfac/4; 4*c/cfac 0 ], ...
                    linear_radius, [ '-' ltraj_color ] );
    
    if include_pause
        pause on;
        disp( 'Press any key to continue' );
        pause
    end
    
% and finally generate solution curves for the nonlinear system
    x_step  = x_max / 6;
    y_step  = y_max / 6;
    x_list  = 0:x_step:x_max;
    y_list  = 0:y_step:y_max;
    ic_list = [ ];
    for j1 = 1:length(x_list)
        ic_list = vertcat( ic_list, [ x_list(j1), y_max ] );
        ic_list = vertcat( ic_list, [ x_list(j1), 0 ] );
    end
    for j2 = 2:length(y_list) - 1
        ic_list = vertcat( ic_list, [ 0, y_list(j2) ] );
        ic_list = vertcat( ic_list, [ x_max, y_list(j2) ] );
    end
    
    for j1 = 1:length(ic_list)
        [t_sol, x_sol] = ode45( @(t,x) [ x(1) + r*x(2) - x(1)*x(2); ...
                                         -c*x(2) + x(1)*x(2) ], ...
                                [ 0 20 ], ic_list(j1,:) );
        plot( figaxes, x_sol(:,1), x_sol(:,2), ['-' ntraj_color], ...
              'LineWidth', 2 );
    end
    
end

% all done
%-------------------------------------------------------------------------------
