%% initialize the system
clear all
% close all
[name,n]=sys_init;
delay_begin = 2;
delay_end = 4;
delay_step = 0.0001;

% start off with an eq. pt.
stst.kind='stst';
stst.parameter=delay_begin; % first parameter value
stst.x=[0 0]';  % approx ic's for eq. pt.
method=df_mthod('stst');
[stst,success]=p_correc(stst,[],[],method.point); % correct if ic's are wrong

% determine stability and plot eigenvalues
method.stability.minimal_real_part=-50; % plot eigenvalues with this large of real part
stst.stability=p_stabil(stst,method.stability);
% figure(1); clf; p_splot(stst) %plot it!

% continue on the steady state branch
branch1 = df_brnch(1,'stst'); % first (and only) parameter (delay)
branch1.parameter.max_bound=[1 delay_end]; % 1 is the parameter pos.
branch1.parameter.max_step = [1 delay_step]; % same as above
branch1.point=stst; % start with the steady state point determined
stst.parameter = delay_begin+delay_step; % modify old stst struct for new delay param
[stst,success]=p_correc(stst,[],[],method.point); % recalculate the stst data

% build up the branch
branch1.point(2)=stst; % next branch point is as calculated
branch1.method.continuation.plot=0; % plot as we're calculating (just stst)
% figure(12); clf;
[branch1,s,f,r]=br_contn(branch1,50000); % the workhorse step, at most 5000 iterates

% stability of branch
branch1 = br_stabl(branch1,0,1); % determine stability of each point on the branch
                                 % note: I edited stst_stabil to disable a
                                 % safeguard against insane corrections.
branch1.point.stability; % usually would print stability of points
branch1.point.parameter; % usually would print delay value of points
%[xm,ym]=df_measr(1,branch1); % get a measure of the eigenvalues along the branch
%figure(4); clf; br_plot(branch1,xm,ym,'b') % plot it
%ym.subfield='l0'; % approximations (l1 is correction)
%
%br_plot(branch1,xm,ym,'c'); % they're nearly indistinguishable
%plot([delay_begin delay_end],[0 0],'-.'); % this is where the axis lies

%% first hopf point
hopf=p_tohopf(branch1.point(10)); % approximate the hopf point to be about here.
method=df_mthod('hopf');
[hopf,success]=p_correc(hopf,1,[],method.point); % correction of hopf pt.
first_hopf=hopf;
hopf.stability=p_stabil(hopf,method.stability); % set the stability of the pt.
figure(3); clf; p_splot(hopf); % plot information on the stability.

% continue non-degenerate periodic solution from first hopf
intervals=20;
degree=3;
[psol,stepcond]=p_topsol(first_hopf,1e-4,degree,intervals); % start with first hopf
method=df_mthod('psol');
[psol,success]=p_correc(psol,1,stepcond,method.point); % correct periodic solution.
branch2=df_brnch(1,'psol');
branch2.parameter.min_bound; % usually would print the minimum bound of the param
branch2.parameter.max_bound=[1 delay_end];
% branch2.parameter.max_step=[1 delay_step]; % default
branch2.parameter.max_step=[1 .01]; % custom
deg_psol=p_topsol(first_hopf,0,degree,intervals); % continue the periodic solution forward
deg_psol.mesh=[]; % save memory by clearing the mesh field
branch2.point=deg_psol;
psol.mesh=[];
branch2.point(2)=psol;
figure(23);
[branch2,s,f,r]=br_contn(branch2,100);

% the following is code to produce profiles of limit cycles
%ll=length(branch2.point);
%figure(10); clf;
%subplot(3,1,1);
%idx = ll-80;
%p_pplot(branch2.point(idx));
%title(['delay = ' num2str(branch2.point(idx).parameter)]);
%xlabel(['t scaled by period = ' num2str(branch2.point(idx).period)]);
%subplot(3,1,2);
%idx = ll-40;
%p_pplot(branch2.point(idx));
%title(['delay = ' num2str(branch2.point(idx).parameter)]);
%xlabel(['t scaled by period = ' num2str(branch2.point(idx).period)]);
%subplot(3,1,3);
%idx = ll-10;
%p_pplot(branch2.point(idx));
%title(['delay = ' num2str(branch2.point(idx).parameter)]);
%xlabel(['t scaled by period = ' num2str(branch2.point(idx).period)]);

%% second hopf point
find_hopf;
hopf = p_tohopf(branch1.point(pt));
method=df_mthod('hopf');
[hopf,success]=p_correc(hopf,1,[],method.point);
second_hopf = hopf;
hopf.stability=p_stabil(hopf,method.stability);

% continue non-degenerate periodic solution from second hopf
intervals=20;
degree=3;
[psol,stepcond]=p_topsol(second_hopf,1e-4,degree,intervals); % start with second hopf
method=df_mthod('psol');
[psol,success]=p_correc(psol,1,stepcond,method.point); % correct periodic solution.
branch3=df_brnch(1,'psol');
branch3.parameter.min_bound; % usually would print the minimum bound of the param
branch3.parameter.max_bound=[1 delay_end];
branch3.parameter.max_step=[1 .01];
deg_psol=p_topsol(second_hopf,0,degree,intervals); % continue the periodic solution forward
deg_psol.mesh=[]; % save memory by clearing the mesh field
branch3.point=deg_psol;
psol.mesh=[];
branch3.point(2)=psol;
figure(23);
[branch3,s,f,r]=br_contn(branch3,100);
