% a demonstration that considers the system
%    x' = [ -r/4 -r/4; 2/5 -2/5 ] x,
% which is a model for the populations of baleen whales (x2)
% and krill (x1), near their scaled equlibrium populations.
%
% if r (= the krill reproductive rate) = 1/3, the eigenvalues and 
% vectors of the coefficient matrix are (1/120)(-29 +- sqrt(119)) and
% [ 19 +- sqrt(119); 48 ].  squinting a little, this suggests taking 
% r = 0.335, for which the eigenvalues are approximately
%    l = -1/4 +- i/11,
% with eigenvectors
%    v = [ 2/5 +- i/4; 1 ]
%
% if r = 4/3, we get a nicer, oscillatory solution; 
%    l = -11/30 +- i sqrt(119)/30 ~ -1/3 +- i/3
%    v = [ 1 +- i sqrt(119); 12 ] ~ [ 1 +- 11 i; 12 ]
%
% by Gavin LaRose, 10/23/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 Whale_Krill_Model()
%-------------------------------------------------------------------------------
% configuration variables

% initial conditions
    x0 = [ 3/4; 3/4 ];
    
% coefficient matrix
    % a_matrix = [ -1/12 -1/12; 2/5 -2/5 ];
    a_matrix = [ -1/3 -1/3; 2/5 -2/5 ];
    
% the maximum time to graph to
    t_max = 20;
    
% should we show the phase plane as well?
    show_phase_plane = 1;
% should there be a pause before showing it?
    include_pause    = 0;
% should we include a sample point on the solutions and phase plane?
    show_point       = 1;
% which point?
    t0 = 3;
    
% colors
    x1_color = 'r';
    x2_color = 'b';
    pp_color = 'k';

    background_grey = [ 0.6 0.6 0.6 ];

%-------------------------------------------------------------------------------
% solution to the problem

% we look at the t range
    t = linspace(0,t_max,1000);

% the eigenvalues and vectors
    [evecs, lambdas] = eig( a_matrix );
% we know these are complex; the real and imaginary parts of the
%     eigenvalues and eigenvectors are
    l_re = real( lambdas(1,1) );
    l_im = imag( lambdas(1,1) );
% the real and imaginary part of an eigenvector are then
    vec1 = real( evecs(1:2,1) );
    vec2 = imag( evecs(1:2,1) );

% we can solve for the constants with
    cvals = [vec1 vec2]^(-1)*x0;

% and our solution is
    % xsol = (cvals(1)*(vec1.*cos(l_im*t) - vec2.*sin(l_im*t)) + ...
    %         cvals(2)*(vec1.*sin(l_im*t) + vec2.*cos(l_im*t))) .* ...
    %        exp(l_re*t);
    x1 = (cvals(1)*(vec1(1)*cos(l_im*t) - vec2(1)*sin(l_im*t)) + ...
          cvals(2)*(vec1(1)*sin(l_im*t) + vec2(1)*cos(l_im*t))) .* ...
         exp(l_re*t);
    x2 = (cvals(1)*(vec1(2)*cos(l_im*t) - vec2(2)*sin(l_im*t)) + ...
          cvals(2)*(vec1(2)*sin(l_im*t) + vec2(2)*cos(l_im*t))) .* ...
         exp(l_re*t);
    
%-------------------------------------------------------------------------------
% graph initialization

% set up the graphics window
    scrsz    = get(0,'ScreenSize');
    scrsz(3) = min(1500, scrsz(3));
    scrsz(4) = min(900, scrsz(4));
    
    if show_phase_plane
        hsc = 1.5;
    else
        hsc = 2.5;
    end
    
    fig1 = figure( 'OuterPosition', ...
                   [50 scrsz(4)/2 scrsz(3)/hsc scrsz(4)/2], ...
                   'Name', 'Model Solution' );
    
    if show_phase_plane
        fig_soln_axes = subplot( 1,2,1, 'Parent', fig1, ...
                                 'Color', background_grey );
        fig_ppln_axes = subplot( 1,2,2, 'Parent', fig1, ...
                                 'Color', background_grey );
        
        hold( fig_ppln_axes );
        xlabel( fig_ppln_axes, 'x1', 'FontSize', 15 );
        ylabel( fig_ppln_axes, 'x2', 'FontSize', 15 );
        title( fig_ppln_axes, 'Phase plane and trajectory', ...
               'FontSize', 15 );
    else
        fig_soln_axes = axes( 'Parent', fig1, 'Color', background_grey );
    end
    
    hold( fig_soln_axes );
    xlabel( fig_soln_axes, 't', 'FontSize', 15 );
    ylabel( fig_soln_axes, 'x1, x2', 'FontSize', 15 );
    title( fig_soln_axes, 'Solution curves', 'FontSize', 15 );

%-------------------------------------------------------------------------------
% plot solution

    sol_plot = plot( fig_soln_axes, ...
                     t, x1, [ x1_color '-' ], ...
                     t, x2, [ x2_color '-' ], 'LineWidth', 2 );
    legend( sol_plot, { 'x1', 'x2' } );
    
    % add axes
    plot( fig_soln_axes, [0 t_max], [0 0], ':k' );
    
    if show_phase_plane
        if include_pause
            pause on;
            disp( 'Press any key to continue' );
            pause;
        end
        
        if show_point
            % the index of the point is
            tind = round((t0*1000)/t_max) + 1;
            plot( fig_soln_axes, ...
                  [t(tind)], [x1(tind)], [ x1_color '.' ], ...
                  [t(tind)], [x2(tind)], [ x2_color '.' ], ...
                  'MarkerSize', 30 );

            plot( fig_ppln_axes, [x1(tind)], [x2(tind)], ...
                  [ pp_color '.' ], 'MarkerSize', 30 );
            
            plot( fig_soln_axes, [0 t0], [x1(tind) x1(tind)], ...
                  [ x1_color '--' ] );
            plot( fig_soln_axes, [0 t0], [x2(tind) x2(tind)], ...
                  [ x2_color '--' ] );
        end
        
        plot( fig_ppln_axes, ...
              x1, x2, [ pp_color '-' ], 'LineWidth', 2 );
        % add an arrow showing the direction of the trajectory
        axes( fig_ppln_axes );
        arrind = round(tind/2) + 1;
        arrow( [ x1(arrind),x2(arrind) ], ...
               [ x1(arrind+1),x2(arrind+1) ] );

        lim = axis;
        plot( fig_ppln_axes, [lim(1) lim(2)], [0 0], ':k' );
        plot( fig_ppln_axes, [0 0], [lim(3) lim(4)], ':k' );
        
        if show_point
            plot( fig_ppln_axes, [x1(tind) x1(tind)], ...
                  [ lim(3) x2(tind) ], [ x1_color '--' ] );
            plot( fig_ppln_axes, [lim(1) x1(tind)], ...
                  [x2(tind) x2(tind)], [ x2_color '--' ] );
        end            
    end
end

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