%
% This script and the file hopftest.m provide a simple
% example of using MATCONT to do one-parameter continuation
% of equilibria and periodic orbits.
%
% The vector field is defined in hopftest.m.
%
% The parameters are mu and omega.
%
% These are the starting parameter values.  mu will be the
% control parameter (i.e. the bifurcation parameter)
%
mu    = -1.5;
omega =  1.0;
p_start = [mu; omega];

%
% eq_start is an equilbrium when mu=-1.5 and omega=1.
% (In fact, (0,0) is an equilibrium for all mu.  This
% makes it easy to test that the code computes the correct
% answer.)
eq_start = [0;0];

%
% ap gives the index of the control parameter
% (mu is the first parameter in the list, so ap=[1].)
%
ap = [1];

%
% Call init_EP_EP to set up MATCONT's global variables
% for continuing an equilibrium solution.  
%
[x0,v0] = init_EP_EP(@hopftest,eq_start,p_start,ap);

%
% Set some continuation options...
%
opt = contset;
opt = contset(opt,'MaxNumPoints',30);
opt = contset(opt,'Singularities',1);

%
% Run the equilibrium continuation
%
[x,v,s,h,f] = cont(@equilibrium,x0,[],opt);

%
% Plot the results
%
cpl(x,v,s,[3 1])
set(gca,'fontsize',14,'fontweight','bold')
xlabel('\mu')
ylabel('x')

%
% Force MATLAB to make the plot, so we have
% something to look at while the limit cycle
% continuation is taking place.
%
drawnow

%
% We should have found a Hopf bifurcation at mu=0.
%
% Print some information about the second point
% recorded in s.  This should be the Hopf point.
%
k = s(2).index;
fprintf('%s at mu=%f,  (%f,%f)\n',s(2).msg,x(3,k),x(1,k),x(2,k));

%
% Now continue the family of periodic orbits from
% the Hopf point.
%
% First initialize MATCONT's global variable with the data
% from the Hopf bifurcation found above by calling init_H_LC
%
% xbifpoint is hold the (x,y) coordinates of the equilibrium
% where the Hopf bifurcation occurs
%
% mubifpoint is the value of mu where the bifurcation occurs
%
xbifpoint = x(1:2,k);
mubifpoint = x(3,k);
ntst = 30;
ncol = 4;
[x0,v0] = init_H_LC(@hopftest,xbifpoint,[mubifpoint; omega],ap,0.001,ntst,ncol);

%
% Set the number of points (i.e. number of periodic orbits)
% to compute.
%
lcnumpoints = 100;
opt = contset;
opt = contset(opt,'MaxNumPoints',lcnumpoints);

%
% Run the limit cycle continuation
%
[xlc,vlc,slc,hlc,flc] = cont(@limitcycle,x0,v0,opt);

%
% Compute the x and y amplitudes of each periodic
% orbit, and then plot the x amplitude as a function
% of the parameter.  (This will be added to the plot
% of the equilibria.)
%
% (There might be another way to get this data.
% I haven't found a MATCONT function that will return
% the k-th periodic orbit in a nice format, but there
% might be one.)
%
xwidth = zeros(lcnumpoints,1);
ywidth = zeros(lcnumpoints,1);
m = 2*(ntst*ncol+1);
for k=1:lcnumpoints,
    xwidth(k) = max(xlc(1:2:(m-1),k))-min(xlc(1:2:(m-1),k));
    ywidth(k) = max(xlc(2:2:m,k))-min(xlc(2:2:m,k));
end

hold on
plot(xlc(m+2,:),xwidth,'g','linewidth',2)
axis([-1.5 1 -0.2 2])