% 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).
%
% a phase portrait for the system near one of the equilibrium solutions
% is graphed.
%
% 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 +/- i 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()
%-------------------------------------------------------------------------------
% configuration variables

% value of c
    c_value = 3;
    
% the equilibrium point to show: 0 = (0,0); 1 = the non-zero equilibrium
    show_point = 0;
    
% should we pause between showing the eigenvalues and vectors (this 
%    is skipped if they are complex), direction field and solution
%    trajectories?  (set to 0 for no; 1 for yes)
    include_pauses = 0;
    
% distance from the equilibria to start trajectories that are attracted
% to the point
    attractor_delta = 0.5;
    
% distance from the equilibria to start trajectories that are moving
% away from the point
    repeller_delta  = 0.01;
    
% length of eigenvector to graph
    arrow_length    = 0.5;
    
% colors
    eigv_color = 'k';  % the color of eigenvectors
    traj_color = 'b';  % the color of solution trajectories
    df_color   = 'k';  % the color of the direction field
    
    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)/2.5 scrsz(3)/2.5], ...
                   'Name', 'Phase Plane' );
    figaxes = axes( 'Parent', fig1, 'Color', background_grey );
    hold(figaxes);
    
% put in some axes
    plot( figaxes, [-.5 .5], [0 0], ':k', ...
          [0 0], [-.5 .5], ':k', 'LineWidth', 2 );

    xlabel( figaxes, 'u', 'FontSize', 15 );
    ylabel( figaxes, 'v', 'FontSize', 15 );

    axis( [ -.5 .5 -.5 .5 ] );
    axis square;

% find the eigenvalues and vectors
    if show_point == 0
        coeff_matrix = [ 1 r; 0 -c ];
    else
        coeff_matrix = [ -9/cfac -cfac/4; 4*c/cfac 0 ];
    end
    [ evecs, lambdas ] = eig( coeff_matrix );
    lam1 = num2str(round(100*lambdas(1,1))/100);
    lam2 = num2str(round(100*lambdas(2,2))/100);
    v11 = num2str(round(100*evecs(1,1))/100);
    v12 = num2str(round(100*evecs(2,1))/100);
    v21 = num2str(round(100*evecs(1,2))/100);
    v22 = num2str(round(100*evecs(2,2))/100);

    if show_point == 0
        title( figaxes, ...
               [ 'Phase portrait near (0,0): \lambda_{1,2} = ', ...
                 lam1 ', ' lam2 '; v_{1,2} = [' v11 '; ' v12 '], [' ...
                 v21 '; ' v22 ']' ], ...
               'FontSize', 15 );
    else
        if imag(lambdas(1,1)) == 0
            title( figaxes, ...
                   [ 'Phase portrait near (' ...
                     num2str(nz_equil(1)) ',' ...
                     num2str(round(100*nz_equil(2))/100) ...
                     '): \lambda_{1,2} = ', ...
                     lam1 ', ' lam2 '; v_{1,2} = [' v11 '; ' v12 '], [' ...
                     v21 '; ' v22 ']' ], ...
                   'FontSize', 15 );
        else
            lam1 = num2str(round(100*real(lambdas(1,1)))/100);
            lam2 = num2str(abs(round(100*imag(lambdas(1,1)))/100));
            
            v11  = num2str(round(100*real(evecs(1,1)))/100);
            v12  = num2str(abs(round(100*imag(evecs(1,1)))/100));
            
            v21  = num2str(round(100*real(evecs(2,1)))/100);
            v22  = num2str(abs(round(100*imag(evecs(2,1)))/100));
            
            if imag(evecs(2,1)) == 0
                title( figaxes, ...
                       [ 'Phase portrait near (' ...
                         num2str(nz_equil(1)) ',' ...
                         num2str(round(100*nz_equil(2))/100) ...
                         '): \lambda = ', lam1 '\pm' lam2 'i; ' ...
                         'v = [' v11 '\pm ' v12 'i; ' v21 ']' ], ...
                       'FontSize', 15 );
            else
                title( figaxes, ...
                       [ 'Phase portrait near (' ...
                         num2str(nz_equil(1)) ',' ...
                         num2str(round(100*nz_equil(2))/100) ...
                         '): \lambda = ', lam1 '\pm' lam2 'i; ' ...
                         'v = [' v11 '; ' v21 '\pm ' v22 'i]' ], ...
                       'FontSize', 15 );
            end
        end
    end
    
%-------------------------------------------------------------------------------
% graphing

% put in the eigenvectors, if eigenvalues are real
    plot_eigenvectors( figaxes, coeff_matrix, arrow_length, eigv_color );
    
    if imag(lambdas(1,1)) == 0 & include_pauses
        pause on;
        disp( 'Press any key to continue' );
        pause
    end
        

% put in the direction field
    plot_directionfield( figaxes, coeff_matrix, 11, df_color );

    if include_pauses
        pause on;
        disp( 'Press any key to continue' );
        pause
    end
    
% put in the solution trajectories
    plot_solutiontraj( figaxes, coeff_matrix, ['-' traj_color] );

end

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