% A demonstration to show the analytical solution to a simplification
% of the system
%    x1' = -(k01 + k21 + k31) x1 + k12 x2 + k13 x3 + I1
%    x2' = k21 x1 - (k02 + k12) x2
%    x3' = k31 x1 - k13 x3
% (the lead modeling problem; c.f. documentation and the demonstration
% Lead_Euler.m).  In this, after substituting values for the constants
% and simplifying to model only the blood and tissues, we can get the 
% system
%    x1' = -(k01 + k21) x1 + k12 x2 + I1
%    x2' = k21 x1 - (k02 + k12) x2
% rescaling a bit, we get the easily solved system
%    x1' = -15 x1 + 4 x2 + I1
%    x2' = 3 x1 - 11 x2,
% If we assume no intake of lead, we can take I1 = 0; then, starting
% at equilibrium, we have, approximately x1(0) = 1300, and x2(0) = 500.
%
% 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 Lead_Model()
%-------------------------------------------------------------------------------
% configuration variables

% initial conditions
    x0 = [ 1300; 500 ];
    
% transfer matrix
    a_matrix = [ -15 4; 3 -11 ];
    
% the maximum time to graph to
    t_max = 1;
    
% 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 = 0.1;
    
% 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 );
% the eigenvectors are then
    vec1 = evecs(1:2,1);
    vec2 = evecs(1:2,2);
% with corresponding eigenvalues
    lam1 = lambdas(1,1);
    lam2 = lambdas(2,2);

% we can solve for the constants with
    cvals = evecs^(-1)*x0;
    
% so that the solution is
    xsol = cvals(1)*vec1*exp(lam1*t) + cvals(2)*vec2*exp(lam2*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, xsol(1,:), [ x1_color '-' ], ...
                     t, xsol(2,:), [ x2_color '-' ], 'LineWidth', 2 );
    legend( sol_plot, { 'x1', 'x2' } );
    
    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 = t0*1000 + 1;
            plot( fig_soln_axes, ...
                  [t(tind)], [xsol(1,tind)], [ x1_color '.' ], ...
                  [t(tind)], [xsol(2,tind)], [ x2_color '.' ], ...
                  'MarkerSize', 30 );

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

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