% predprey_bif.m - bifurcation plots for predator prey analysis
% RMM, 8 Sep 06

aminit;					% initialize macros

% Load Mathematica data from file
load predprey.dat

% Store variables in easier places
a = predprey(:,1);
Heq = predprey(:,2);
stab = predprey(:,3);
Hmn = predprey(:,4);
Hmx = predprey(:,5);
cmin = predprey(:,6);
cmax = predprey(:,7);

% Tweak the cmax values
starti = find(cmax >= 200, 1, 'first');
cmax(starti:length(cmax)) = 200;

% Generate the plot for the stability boundary
figure(1); clf; subplot(221);  hold on;
patch([a; a(length(a)); a(starti); a(1)], [cmin; 200; 200; cmax(1)], ...
  [0.9 0.9 0.9]);
plot([4 4], [0 200], 'k', 'LineWidth', 1);	% fix box border
amaxis([1.35, 4, 0, 200]); box on;
xlabel('{\ita}');
ylabel('{\itc}', 'Rotation', 0);
amtext(1.45, 170, 'Unstable');
amtext(2.25, 100, 'Stable');
amtext(3, 25, 'Unstable');

amprint('predprey-parm.eps');

% Go through the data and sort out stable versus unstable
segi = 1; ind = 1; signx = sign(stab(1));
seglen(1) = 0;
for i = 1:length(a)
  if (stab(i)*signx < 0) 
    % Stability changed
    segi = segi + 1;
    ind = 1;
    signx = sign(stab(i));
    seglen(segi) = 0;
  end

  % Store the data
  aseg(segi, ind) = a(i);
  Heqseg(segi, ind) = Heq(i);
  stabseg(segi, ind) = stab(i);

  % Increment pointers
  seglen(segi) = seglen(segi) + 1;
  ind = ind + 1;
end
numseg = segi;

% Plot the bifurcation diagram
figure(2); clf; subplot(221);
for segi = 1:numseg
  if (stabseg(segi, 1) < 0) 
    plot(aseg(segi, 1:seglen(segi)), Heqseg(segi, 1:seglen(segi)), 'b-', ...
    'Linewidth', AM_thick_linewidth);
  else
    plot(aseg(segi, 1:seglen(segi)), Heqseg(segi, 1:seglen(segi)), 'r--', ...
    'Linewidth', AM_thick_linewidth);
  end
  hold on;
end

% Plot the limit cycle bounds
plot(a, Hmn, 'k-.', a, Hmx, 'k-.', 'Linewidth', AM_data_linewidth);
amaxis([1.35 8 0 150]);
xlabel('{\ita}');
ylabel('{\itH}', 'Rotation', 0);

amprint('predprey-bif.eps');
