% A demonstration to show the phase plane and solution to the system
%    x1' = I - (k1 + k2) x1
%    x2' = k1 x1 - k3 x2.
% Because we need the system to be autonomous to plot a phase plane, we 
% consider I = 0 and take a non-zero initial condition for the system.
% 
% Alternately, we consider another function, as defined in the function
% rhs(t,x).)  A configuration variable determines whether the solution 
% and phase plane are subplots of a single graph window or to different
% windows.  A second configuration variable allows an easy conversion
% of the rhs function to give an oscillatory solution.  This second
% system is 
%    x1' = -x1 - 2 x2
%    x2' = 6 x1
% initial conditions in this case need to be non-zero; in the case of
% the drug modeling problem they may be zero.
%
% by Justin Dunmyre
% revised by Gavin LaRose, 10/12/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 Phase_Plane()
%-------------------------------------------------------------------------------
% configuration variables

    k1 = 0.128;
    k2 = 0.383;
    k3 = 0.0385;
    I0 = 25;
    
    % the initial conditions: the amount of the drug in each
    %    compartment initially
    x0 = [ I0; 0 ];
    
    % should we graph the solution and phase plane in different
    %    windows?  (set to 0 for no, 1 for yes)
    use_two_windows = 0;
    
    % should we use the alternate, oscillatory, system solution?
    use_alt_solution = 0;
    
    % in this case we reset the initial conditions
    if use_alt_solution == 1
        x0 = [ 4; 0 ];
    end
    
    % the range of times for which to show the solution
    if use_alt_solution == 1
        t_range = [ 0 20 ];
    else
        t_range = [ 0 120 ];
    end
    
    % animation speed
    if use_alt_solution == 1
	speed = 1;
    else
	speed = 2;
    end
    
    % colors: x1, x2, the phase plane trajectory, and the direction field
    x1_color = 'b';
    x2_color = 'r';
    pp_color = 'g';
    df_color = 'k';
    
    background_grey = [ 0.6 0.6 0.6 ];
    
    % number of vectors shown in each direction in the direction field
    num_vec = 11;
    
%-------------------------------------------------------------------------------
% system definition: a function for the right-hand side of the system

    function out = rhs( t, x )
        if use_alt_solution == 1
            out = [ -x(1) - 2*x(2); 6*x(1) ];
        else
            out = [ -(k1+k2)*x(1); k1*x(1) - k3*x(2) ];
        end
    end 
    
%-------------------------------------------------------------------------------
% initialization
    
% set up the graphics 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(4)/2], ...
                   'Name', 'Model Solution' );
    fig1_axes = axes( 'Parent', fig1, 'Color', background_grey );
    % this ensures we have a toolbar
    set( fig1, 'Toolbar', 'figure' );
    
% if we're setting two plot windows, we generate a second figure
%    window; otherwise, we split this into two subplots
    if use_two_windows == 1
        fig2 = figure( 'OuterPosition', ...
                       [100 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2], ...
                       'Name', 'Phase Plane' );
        set( fig2, 'Toolbar', 'figure' );
        
        fig_ppln_axes = axes( 'Parent', fig2, 'Color', background_grey );
        fig_soln_axes = fig1_axes;
    else
        fig_soln_axes = subplot( 1,2,1, 'Parent', fig1, ...
                                 'Color', background_grey );
        fig_ppln_axes = subplot( 1,2,2, 'Parent', fig1, ...
                                 'Color', background_grey );
    end
    
    hold( fig_soln_axes );
    hold( fig_ppln_axes );
    xlabel( fig_soln_axes, 't', 'FontSize', 15 );
    ylabel( fig_soln_axes, 'x1, x2', 'FontSize', 15 );
    xlabel( fig_ppln_axes, 'x1', 'FontSize', 15 );
    ylabel( fig_ppln_axes, 'x2', 'FontSize', 15 );
    title( fig_soln_axes, 'Solution curves', 'FontSize', 15 );
    title( fig_ppln_axes, 'Phase plane, direction field and trajectory', ...
           'FontSize', 15 );
    
%-------------------------------------------------------------------------------
% solution and direction field generation and initial plot

    [t_sol, x_sol] = ode45( @rhs, linspace(t_range(1),t_range(2),100), x0 );

    % the extreme values for x1 and x2 are 
    x1_max = max( x_sol(:,1) );
    x2_max = max( x_sol(:,2) );
    x1_min = min( x_sol(:,1) );
    x2_min = min( x_sol(:,2) );
    x_max  = round(max( [ x1_max x2_max ] ));
    x_min  = round(max( [ x1_min x2_min ] ));
    
    % set the size of the axes
    axis( fig_soln_axes, [ t_range(1), t_range(2), x_min-1, x_max+1 ] );
    axis( fig_ppln_axes, [ x1_min-1, x1_max+1, x2_min-1, x2_max+1 ] );
    
    % set up x and y vectors for the direction field
    x1_vec = linspace((x1_min-1),(x1_max + 1), num_vec);
    x2_vec = linspace((x2_min-1),(x2_max + 1), num_vec);
    
    % then we define an x and y grid and the corresponding derivatives
    x1_grid = [];
    x2_grid = [];
    x1p_grid = [];
    x2p_grid = []; % zeros(num_vec, num_vec);
    
    for j1=1:num_vec
        % x1_grid( j1,: ) = x1_vec;
        % x2_grid( :,j1 ) = x2_vec;
        for j2=1:num_vec
            
            x1_grid = [ x1_grid, x1_vec(j1) ];
            x2_grid = [ x2_grid, x2_vec(j2) ];

            der = rhs( 0, [ x1_vec(j1); x2_vec(j2) ] );
            % x1p_grid(j1,j2) = der(1);
            % x2p_grid(j1,j2) = der(2);
            x1p_grid = [ x1p_grid der(1) ];
            x2p_grid = [ x2p_grid der(2) ];
        end
    end

    % [ x_vec, y_vec ] = meshgrid( (-x_max_abs-1):dx:(x_max_abs + 1) );

    % x_vec(1,:) row vec of x values
    % y_vec(:,1) col vec of y values
    
    % to normalize arrow lengths, the lengths of these are
    dx_len = sqrt(x1p_grid.^2 + x2p_grid.^2);
    
    % then the direction field is
    quiver( fig_ppln_axes, ...
            x1_grid, x2_grid, x1p_grid./dx_len, x2p_grid./dx_len, ...
            'ShowArrowHead', 'on', 'AutoScaleFactor', .3, 'Color', df_color );
    
    % plot the solutions with dotted lines
    plot( fig_soln_axes, ...
          t_sol, x_sol(:,1), [ x1_color ':' ], ...
          t_sol, x_sol(:,2), [ x2_color ':' ], ...
          'LineWidth', 2 );
    
    % put a dot on the initial condition
    xplots = plot( fig_soln_axes, ...
                   0, x0(1), [ x1_color '.' ], ...
                   0, x0(2), [ x2_color '.' ], 'MarkerSize', 20 );
    legend( xplots, { 'x1(t)', 'x2(t)' } );

    pplot = plot( fig_ppln_axes, x0(1), x0(2), [ pp_color '.' ], ...
                  'MarkerSize', 20 );

    getframe();
  
%-------------------------------------------------------------------------------
% set up the animation

    % initialize containers for the graphs
    phase_dot_gr  = zeros(length(t_sol), 1);
    phase_traj_gr = zeros(length(t_sol), 1);
    
    x1_dot_gr  = zeros(length(t_sol), 1);
    x2_dot_gr  = zeros(length(t_sol), 1);
    x1_soln_gr = zeros(length(t_sol), 1);
    x1_soln_gr = zeros(length(t_sol), 1);

    for j=1:length(t_sol)
        phase_dot_gr(j)  = plot( fig_ppln_axes, ...
                                 x_sol(j,1), x_sol(j,2), ...
                                 [ pp_color '.' ], 'MarkerSize', 20, ...
                                 'Visible', 'off' );
        phase_traj_gr(j) = plot( fig_ppln_axes, ...
                                 x_sol(1:j,1), x_sol(1:j,2), ...
                                 [ pp_color '-' ], 'LineWidth', 2, ...
                                 'Visible', 'off' );
        x1_dot_gr(j)  = plot( fig_soln_axes, ...
                              t_sol(j), x_sol(j,1), [ x1_color '.' ], ...
                              'MarkerSize', 20, 'Visible', 'off' );
        x2_dot_gr(j)  = plot( fig_soln_axes, ...
                              t_sol(j), x_sol(j,2), [ x2_color '.' ], ...
                              'MarkerSize', 20, 'Visible', 'off' );
        x1_soln_gr(j) = plot( fig_soln_axes, ...
                              t_sol(1:j), x_sol(1:j,1), [ x1_color '-' ], ...
                              'LineWidth', 2, 'Visible', 'off' );
        x2_soln_gr(j) = plot( fig_soln_axes, ...
                              t_sol(1:j), x_sol(1:j,2), [ x2_color '-' ], ...
                              'LineWidth', 2, 'Visible', 'off');
    end

    disp( 'Press any key to continue' );
    pause on;
    pause;
    
    set( xplots, 'Visible', 'off' );
    set( pplot,  'Visible', 'off' );
    
    % then do the animation graphing
    for j=(1+speed):speed:length(t_sol)
        set( phase_dot_gr(j-speed), 'Visible', 'off' );
        set( x1_dot_gr(j-speed), 'Visible', 'off' );
        set( x2_dot_gr(j-speed), 'Visible', 'off' );
        
        set( phase_dot_gr(j), 'Visible', 'on' );
        set( x1_dot_gr(j), 'Visible', 'on' );
        set( x2_dot_gr(j), 'Visible', 'on' );
        
        set( phase_traj_gr(j), 'Visible', 'on' );
        set( x1_soln_gr(j), 'Visible', 'on' );
        set( x2_soln_gr(j), 'Visible', 'on' );
        
        getframe();
    end
end
% all done
%-------------------------------------------------------------------------------
