% Modeling the phenomenon of spruce budworm epidemics
% 
% This demo graphs the right-hand-side of the differential equation
% modeling the budworm population, 
%    p' = r p ( 1 - p/q ) - p^2/(1 + p^2),
% with r = 0.5 and different values of q, and then shows the
% corresponding phase diagrams.
% 
% (c)2013 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_Bifurcation_Sequence()

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

% value of q to use
%    these are the values used in the other demos for this section
%    changing these will require changing the subplot definitions and
%    the phase diagram plots, below.
%   interestingQ=[10 9 8 7.3878 7 6];
    interestingQ=[8 7.3878 7];
    
    % do we pause after showing the first right-hand side?  (1 = yes; 0 = no)
    pause_first = 1;
% do we pause after the right-hand sides and before showing each phase
%    diagram?  (1 = yes; 0 = no)
    pause_diagrams = 1;

% 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'};
    colorvec={'w','m','y'};

% the value of r
    r = 0.5;
    
% more color
    background_grey = [ 0.6 0.6 0.6 ];

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

%--------------------
% set up the figure and plots
    % figh=figure('Position',[scrsz(3)/4 scrsz(4)/4 scrsz(3)/4 scrsz(4)/2]);
    figh = figure();
    subp1=subplot(6,1,[1 3], 'Parent', figh, 'Color', background_grey );
    axis( [ 0 7 -0.25 0.25 ] );
    
    subp2=subplot(6,1,4, 'Parent', figh, 'Color', background_grey);
    axis( [ 0 7 -0.25 0.25 ] );
    subp3=subplot(6,1,5, 'Parent', figh, 'Color', background_grey);
    axis( [ 0 7 -0.25 0.25 ] );
    subp4=subplot(6,1,6, 'Parent', figh, 'Color', background_grey);
    axis( [ 0 7 -0.25 0.25 ] );
    hold(subp1);  hold(subp2);  hold(subp3);  hold(subp4);
    
    title( [ 'Right-hand side of P'' = ' num2str(r) 'P(1 - P/Q) - ' ...
             'P^2/(1 + P^2)' ], 'Parent', subp1, 'FontSize', 16 );
    % xlabel('P', 'Parent', subp1, 'FontSize', 14);
    ylabel('RHS(P)', 'Parent', subp1, 'FontSize', 14);

    title( 'Phase diagrams', 'Parent', subp2, 'FontSize', 16 );

    q = interestingQ(1);
    pl = zeros(1, length(interestingQ));
    leg = cell(1, length(interestingQ));

    % plot the right-hand side in the top plot
    pvals = 0:.01:7;
    pl(1) = plot( pvals, r*pvals.*(1 - pvals/q) - pvals.^2./(1 + pvals.^2), ...
                  [colorvec{1} '-'], 'LineWidth', 2, 'Parent', subp1 );
    leg{1} = ['q = ' num2str(interestingQ(1)) ];
    lh = legend( subp1, leg{1} );
    plot( pvals, zeros(1,length(pvals)), '-k', 'Parent', subp1 );

    if pause_first
        pause on;
        disp( 'Press enter for next graphs' );
        pause;
    end

    for i = 2:length(interestingQ)
        q = interestingQ(i);
        pl(i) = plot( pvals, r*pvals.*(1 - pvals/q) - ...
                      pvals.^2./(1 + pvals.^2), ...
                      [colorvec{i} '-'], 'LineWidth', 2, 'Parent', subp1 );
        leg{i} = [ 'q = ' num2str(interestingQ(i)) ];
    end
    delete(lh);
    lh = legend( pl, leg );

    % then we graph the phase diagrams in the lower plots
    if pause_diagrams
        disp( 'Press enter for phase diagrams' );
        pause;
    end
    
    % first phase diagram
    c = colorvec{1};
    q = interestingQ(1);
    plot( pvals, zeros(1,length(pvals)), ['-' c], ...
          .25, 0, ['>' c], 1.5, 0, ['<' c], 3.5, 0, ['>' c], ...
          7, 0, ['<' c], ...
          'MarkerSize', 15, 'MarkerFaceColor', c, 'Parent', subp2 );
    plot( 0, 0, ['.' c], 0.66, 0, ['.' c], 2.53, 0, ['.' c], ...
          4.81, 0, ['.' c], ...
          'MarkerSize', 30, 'MarkerFaceColor', c, 'Parent', subp2 );

    ylabel( ['q = ', num2str(q)], 'Parent', subp2, 'FontSize', 14 );

    if pause_diagrams
        disp( 'Press enter for next phase diagram' );
        pause;
    end

    % second phase diagram
    c = colorvec{2};
    q = interestingQ(2);
    plot( pvals, zeros(1,length(pvals)), ['-' c], ...
          .25, 0, ['>' c], 2, 0, ['<' c], 6, 0, ['<' c], ...
          'MarkerSize', 15, 'MarkerFaceColor', c, 'Parent', subp3 );
    plot( 0, 0, ['.' c], 0.64, 0, ['.' c], 3.37, 0, ['.' c], ...
          'MarkerSize', 30, 'MarkerFaceColor', c, 'Parent', subp3 );

    ylabel( ['q = ', num2str(q)], 'Parent', subp3, 'FontSize', 14 );

    % second phase diagram
    c = colorvec{3};
    q = interestingQ(3);
    plot( pvals, zeros(1,length(pvals)), ['-' c], ...
          .25, 0, ['>' c], 4, 0, ['<' c], ...
          'MarkerSize', 15, 'MarkerFaceColor', c, 'Parent', subp4 );
    plot( 0, 0, ['.' c], 0.64, 0, ['.' c], ...
          'MarkerSize', 30, 'MarkerFaceColor', c, 'Parent', subp4 );

    ylabel( ['q = ', num2str(q)], 'Parent', subp4, 'FontSize', 14 );

end

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