% A demonstration to show the solution to the system
%    x1' = I - (k1 + k2) x1
%    x2' = k1 x1 - k3 x2,
% where I = (I0/h) (1 - step(t-h)).  Several configuration variables 
% determine how the solution is graphed---whether the input and
% response appear on different axes, and whether we should repeat the
% input every r hours.
%
% 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 Drug_Model()
%-------------------------------------------------------------------------------
% configuration variables

    k1 = 0.128;
    k2 = 0.383;
    k3 = 0.0385;
    I0 = 25;
    h  = 3;
    
    % the initial conditions: the amount of the drug in each
    %    compartment initially
    x0 = [ 0; 0 ];
    
    % should we repeat the input? (set to 0 for no, 1 for yes)
    repeat_input = 1;
    % how often?  every repeat_time hours:
    repeat_time  = 6;
    
    % should we graph the input and response on the same axes, or
    %    different?  (set to 0 for no, 1 for yes)
    use_same_axes = 1;
    
    % the range of times for which to show the solution
    t_range = [ 0 48 ];
    
    % colors
    x1_color  = 'b';
    x2_color  = 'r';
    inp_color = 'k';
    
    background_grey = [ 0.6 0.6 0.6 ];
    
%-------------------------------------------------------------------------------
% initialization
    
% set up the graphics window
    scrsz    = get(0,'ScreenSize');
    scrsz(3) = min(1500, scrsz(3));
    scrsz(4) = min(900, scrsz(4));

    if use_same_axes == 0
        h_scale = 1.5;
    else
        h_scale = 2;
    end
    figh = figure( 'OuterPosition', ...
                   [50 scrsz(4)/2 scrsz(3)/2 scrsz(4)/h_scale], ...
                   'Name', 'Model Solution' );
    fig_axes = axes( 'Parent', figh, 'Color', background_grey );
    % this ensures we have a toolbar
    set( figh, 'Toolbar', 'figure' );

% if we're setting up two plots, do that; otherwise we're ok with one
    if use_same_axes == 0
        soln_axes  = subplot( 2,1,1, 'Parent', figh, 'Color', background_grey );
        input_axes = subplot( 2,1,2, 'Parent', figh, 'Color', background_grey );
        set( soln_axes, 'xlim', t_range );
        set( input_axes, 'xlim', t_range );
        
        % set(soln_plot, 'Xtick', t_range(1):1:t_range(2));       
        % set(input_plot, 'Xtick', t_range(1):1:t_range(2));       
        
        hold( soln_axes );
        hold( input_axes );       

        xlabel( soln_axes, 't (h)', 'FontSize', 15 );
        xlabel( input_axes, 't (h)', 'FontSize', 15 );
        ylabel( soln_axes, 'drug (mg)', 'FontSize', 15 );
        ylabel( input_axes, 'drug (mg)', 'FontSize', 15 );
    
        title( soln_axes, 'Drug response', 'FontSize', 15 );
        title( input_axes, 'Input function', 'FontSize', 15 );
    else
        soln_axes  = fig_axes;
        input_axes = fig_axes;
        
        set( fig_axes, 'xlim', t_range );

        hold on;
        
        xlabel( 't (h)', 'FontSize', 15 );
        ylabel( 'drug (mg)', 'FontSize', 15 );
        title( 'Input and drug response', 'FontSize', 15 );
    end        
    
% set up a function for the input
    function out = dose(t)
        out = (I0/h)*(t < h | ( repeat_input & rem(t,repeat_time) < h ));
    end
    
% set up a function for the right-hand side of the system
    function out = rhs( t, x )
        out = [ dose(t) - (k1+k2)*x(1); k1*x(1) - k3*x(2) ];
    end
    
% then the solution to the system is
    [ t_sol, x_sol ] = ode45( @rhs, t_range, x0 );
    
%-------------------------------------------------------------------------------
% plot the solution

    xplots = plot( soln_axes, ...
                   t_sol, x_sol(:,1), [ x1_color '-' ], ...
                   t_sol, x_sol(:,2), [ x2_color '-' ], ...
                   'LineWidth', 2 );
    iplot  = plot( input_axes, ...
                   t_sol, dose(t_sol), [ inp_color '-' ], ...
                   'LineWidth', 2 );
    
    if use_same_axes == 0
        legend( xplots, { 'x1(t)', 'x2(t)' } );
        legend( iplot, { 'I(t)' } );
    else
        legend( 'x1(t)', 'x2(t)', 'I(t)' );
    end

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