% Modeling the phenomenon of spruce budworm epidemics

% This demo is just the last portion of the Budworm_Model
%    demonstration, showing the animated phase plane only.
% by Justin Dunmyre
% modified by Gavin LaRose, 8/10/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 budworm()

%----------------------------------------------------------------------
% configuration variables

% value of q to use
%    these are the values used in the other demos for this section
% interestingQ=[10 9 8 7.3878 7 6];
    qval = 10;

% the corresponding colors we use for these solutions
%    these are the values used in the other demos
    colorvec={'w','m','y','c','g','b','k'};
    col = 'w';

% an initial condition
    ic = 3;

% the value of r
    r = 0.5;

%----------------------------------------
% the rest of this should not need modification

%--------------------
% the right-hand side of the differential equation
    function uprime=budwormODE(t,u)
	uprime=r.*u.*(1-u/q)-((u.^2)./(1+u.^2));
    end

%--------------------
% set up the figure and plots
    % figh=figure('Position',[scrsz(3)/4 scrsz(4)/4 scrsz(3)/4 scrsz(4)/2]);
    figh = figure();
    subplot1=subplot(2,1,1, 'Parent', figh, 'Color',[0.5, 0.5, 0.5]);
    subplot2=subplot(2,1,2, 'Parent', figh, 'Color',[0.5, 0.5, 0.5]);
    hold(subplot1);
    hold(subplot2);

    q = qval;
    [t,u]=ode45(@budwormODE,[0 200], ic);

    % the solution curve is graphed in the top figure
    plot( [t(1), t(1)+.01], [u(1), u(1)+.01*((u(2)-u(1))/(t(2)-t(1)))],...
        [col '-'], 'LineWidth', 2, 'Parent',subplot1);
    legend(subplot1, {'Solution to ODE'}, 'location','SE')

    p1 = plot( t(1), u(1), [col '.'],...
	    'MarkerSize', 30, 'Parent',subplot1);
    plot(t,u,[col ':'], 'LineWidth', 2, 'Parent', subplot1) 
    xlabel('t', 'Parent', subplot1);
    ylabel('u', 'Parent', subplot1);

    % then we graph the phase plane in the lower figure
    x=linspace(0,qval+1,10000);
    plot(x,budwormODE(1,x),'k-','LineWidth',2,'Parent',subplot2);
	plot( [u(1), u(1)+.1], [0, 0], [col '-'],...
	      'LineWidth', 2, 'Parent',subplot2);
    legend(subplot2,{'RHS of ODE','Solution to ODE'})

    p2 = plot( u(1), 0, [col '.'], 'MarkerSize', 30, 'Parent',subplot2);

    % put in the axis as a dotted line
    plot(x,zeros(1,length(x)),'k:','LineWidth',2,'Parent',subplot2);
    
    xlabel('u');
        
    legend(subplot2,{'RHS of ODE','Solution to ODE'})
    axis([0 max(x) -.375 .375]);
    annotation(figh,'textbox', [0.54 0.70 0.17 0.09],...
        'String',{['q=' num2str(q)]},'FontSize',20,...
        'FitBoxToText','on');

    PhaseLineArrowsu=[.25, 1.3, 3.5, 9];
    for j=1:length(PhaseLineArrowsu)
        myu=PhaseLineArrowsu(j);
        if(budwormODE(1,myu)>0)
            PhaseArrow(j)=plot(myu,0,'k>','MarkerSize',15,'MarkerFaceColor','k');
        else
            PhaseArrow(j)=plot(myu,0,'k<','MarkerSize',15,'MarkerFaceColor','k');
        end
    end
    % add the equilibria
    eqn = [ num2str(r), '*x*(1 - x/', num2str(qval), ') - ((x^2)/(1 + x^2))' ];
    evals = solve( eqn );
    for j=1:length(evals)
        plot( evals(j), 0, 'k.', 'MarkerSize', 30, 'MarkerFaceColor', ...
	    'k', 'Parent', subplot2 );
    end

    pause

    % delete the points showing the initial conditions
    delete(p1);
    delete(p2);

    % do the animation by plotting an increasing percentage of
    %   the solution vector appropriately
    for j=1:length(t)
        tvec=t(1:j);
        uvec=u(1:j);
        solutionDotX=t(j);
        solutionDotY=u(j);

        phaseVecX=u(1:j);
        phaseVecY=zeros(1,j);
        phaseDotX=u(j);
        phaseDotY=0;

        plot( tvec, uvec, [col '-'], 'LineWidth', 2, 'Parent', subplot1 );
        soldot = plot( solutionDotX, solutionDotY, [col '.'], ...
            'MarkerSize',30, 'Parent',subplot1 );

        plot( phaseVecX, phaseVecY, [col '-'], 'LineWidth', 2, ...
            'Parent', subplot2);
        phasedot = plot( phaseDotX, phaseDotY, [col '.'], 'MarkerSize', 30,...
            'Parent', subplot2 );
        
        getframe();
        delete(soldot);
        delete(phasedot);
    end
end

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