%% This script reproduces the figures as in the tutorial for the GUI
% We assume the system has been created, and now we need the handle to the
% ODEs for simulations
init;fun=feval('JR_slow');
ODE =fun{2};
%% Simulations for two different parameter values
% setting the parameter
% II=20,A=3.25,a=100,Bs=8.8,bs=20,Bf=44,bf=100,e0=2.5,r=0.56,v0=6,theta=0.5,CC=220;
disp('Simulations');
p=[20;3.25;100;8.8;20;44;100;2.5;0.56;6;0.5;220]; 
% par_II,par_A,par_a,par_Bs,par_bs,par_Bf,par_bf,par_e0,par_r,par_v0,par_theta,par_CC
pp=num2cell(p);
[t,x] = ode45(ODE,[0 1],zeros(8,1),[],pp{:});
figure(20);
subplot(1,2,1);plot(t,x(:,1));
xlabel('Time [s]','Fontsize',14);ylabel('y_{0}','Fontsize',14);
subplot(1,2,2);plot(x(:,1),x(:,3));
xlabel('y_{0}','Fontsize',14);ylabel('y_{1}','Fontsize',14);
p(1)=135;
pp=num2cell(p);
[tp,xp] = ode45(ODE,[0 1],zeros(8,1),[],pp{:});
figure(21);
subplot(1,2,1);plot(tp,xp(:,1));
xlabel('Time [s]','Fontsize',14);ylabel('y_{0}','Fontsize',14);
subplot(1,2,2);plot(xp(:,1),xp(:,3));
xlabel('y_{0}','Fontsize',14);ylabel('y_{1}','Fontsize',14);

%% Continuation of Equilibria in one Parameter
disp('Continuation of Equilibria');
[x0,v0]=init_EP_EP(@JR_slow,x(end,:)',p, 1);
opt=contset;
opt=contset(opt,'MaxNumPoints',600);
opt=contset(opt,'MaxStepsize',5);
opt=contset(opt,'Singularities',1);
opt=contset(opt,'Backward',0);
[x1,v1,s1,h1,f1]=cont(@equilibrium,x0,[],opt);
opt=contset(opt,'Backward',1);
[x2,v2,s2,h2,f2]=cont(@equilibrium,x0,[],opt);
figure(22);clf(22);hold on;
cpl(x1,v1,s1,[9 1]);
cpl(x2,v2,s2,[9 1]);
xlabel('II');ylabel('y_0');xlim([-250 300]);

%% Continuation of Codimension 1 Bifurcations; Limit Point and Hopf curves
%First the Limit Point curve
disp('Continuation of LP and Hopf');
ind1=s1(2).index;ap=[1 12];
x0=x1(1:end-1,ind1);p1=p;p1(1)=x1(end,ind1);
[x0,v0]=init_LP_LP(@JR_slow,x0,p1,ap);
opt=contset(opt,'MaxStepsize',1); % Otherwise we miss the BT points
opt=contset(opt,'MaxNumPoints',800);%To reach the end of the plotting region
opt=contset(opt,'Backward',0);
[x3,v3,s3,h3,f3]=cont(@limitpoint,x0,[],opt);
opt=contset(opt,'Backward',1);
[x4,v4,s4,h4,f4]=cont(@limitpoint,x0,[],opt);
figure(23);clf(23);hold on;
cpl(x3,v3,s3,[9 10]);
cpl(x4,v4,s4,[9 10]);
xlabel('II');ylabel('CC');xlim([-400 300]);ylim([0 400]);
% Now the two Hopf bifurcations
opt=contset(opt,'MaxNumPoints',1800);
ind2=s1(6).index;ap=[1 12]; % First Hopf
x0=x1(1:end-1,ind2);p2=p;p2(1)=x1(end,ind2);
[x0,v0]=init_H_H(@JR_slow,x0,p1,ap);
opt=contset(opt,'Backward',0);
[x5,v5,s5,h5,f5]=cont(@hopf,x0,[],opt);
opt=contset(opt,'Backward',1);
[x6,v6,s6,h6,f6]=cont(@hopf,x0,[],opt);
ind2=s1(7).index;ap=[1 12]; %Second Hopf
x0=x1(1:end-1,ind2);p2=p;p2(1)=x1(end,ind2);
[x0,v0]=init_H_H(@JR_slow,x0,p1,ap);
opt=contset(opt,'Backward',0);
[x7,v7,s7,h7,f7]=cont(@hopf,x0,[],opt);
opt=contset(opt,'Backward',1);
[x8,v8,s8,h8,f8]=cont(@hopf,x0,[],opt);

cpl(x5,v5,s5,[9 10]);cpl(x6,v6,s6,[9 10]);
cpl(x7,v7,s7,[9 10]);cpl(x8,v8,s8,[9 10]);

%% Continuation of Limit Cycles
%Part 1; Starting from a Hopf bifurcation
disp('Continuation of Limit Cycles')
ind2=s1(7).index;ap=[1]; % First Hopf
x0=x1(1:end-1,ind2);p2=p;p2(1)=x1(end,ind2);
[x0,v0]=init_H_LC(@JR_slow,x0,p2,ap,1e-2,40,4);
opt=contset(opt,'MaxNumPoints',300);
opt=contset(opt,'MaxStepsize',100);
opt=contset(opt,'Multipliers',1);
opt=contset(opt,'Adapt',1);
opt=contset(opt,'IgnoreSingularity',[1 4]);
opt=contset(opt,'Backward',0);
[xlc1,vlc1,slc1,hlc1,flc1]=cont(@limitcycle,x0,v0,opt);
figure(24);hold on;
plotcycle(xlc1,vlc1,slc1,[size(xlc1,1) 1 3]);
xlabel('II');ylabel('y_0');zlabel('y_1');
title('A one-parameter family of limit cycles');

%% Part 2; Starting from a simulation
%Take earlier simulation and just "redo" it to be sure.
p(1)=135;pp=num2cell(p);ap=[1];ntst=40;ncol=4;tolerance=1e-2;
[ ~,xp]=ode45(ODE,[0 1],zeros(8,1),[],pp{:});
%Now simulate for a time span lasting longer than one period
[tp,xp]=ode45(ODE,[0 0.4],xp(end,:),[],pp{:});
[x0,v0]=initOrbLC(@JR_slow, tp, xp, p, ap, ntst, ncol, tolerance);
opt=contset;
opt=contset(opt,'Multipliers',1);
opt=contset(opt,'MaxNumPoints',500);
opt=contset(opt,'MaxStepsize',100);
opt=contset(opt,'Singularities',1);
opt=contset(opt,'IgnoreSingularity',[1 4]); %Ignore everything except LPC and PD
opt=contset(opt,'Backward',0);
[xlc2,vlc2,slc2,hlc2,flc2]=cont(@limitcycle,x0,v0,opt) ;
opt=contset (opt,'Backward',1);
[x0,v0]=initOrbLC(@JR_slow, tp, xp, p, ap, ntst, ncol, tolerance);
[xlc3,vlc3,slc3,hlc3,flc3]=cont(@limitcycle,x0,v0,opt) ;
figure(25);plotcycle(xlc2,vlc2,slc2,[size(xlc2,1) 1 3]);
% plotcycle(xlc3,vlc3,slc3,[size(xlc3,1) 1 3]);
xlim([100 200]);
%% Continuation of a Limit Point of Cycles Bifurcation in two parameters.
% We've taken really large steps, so if just start with the continuation of
% the LPC, the initial data does not converge. So we start from the
% detected LPC to improve its detection.
disp('Continuation of Limit Point of Cycles')
opt=contset(opt,'MaxNumPoints',10);
opt=contset(opt,'Singularities',1);
opt=contset(opt,'MaxStepsize',.1);
opt=contset(opt,'Backward',1);
[x0,v0]=init_LC_LC(@JR_slow,xlc2,vlc2,slc2(5),p,[1],40,4);
[xlc4,vlc4,slc4,hlc4,flc4]=cont(@limitcycle,x0,v0,opt) ;
% Now that point has been refined, and we proceed with the LPC continuation.
[x0,v0]=init_LPC_LPC(@JR_slow,xlc4,slc4(2),[1 12],40,4);
opt=contset(opt,'MaxNumPoints',1000);
opt=contset(opt,'Singularities',0);
opt=contset(opt,'Backward',0);
opt=contset(opt,'MaxStepsize',15);
[x9,v9,s9,h9,f9]=cont(@limitpointcycle,x0,[],opt);
figure(23);plot(x9(end-1,:),x9(end,:));
% May be nice to observe this curve ends up near a Generalized Hopf (GH)
% which is expected.
% This LPC continuation is slow, so for now it is commented.
% opt=contset (opt,'Backward',1);
% [x10,v10,s10,h10,f10]=cont(@limitpointcycle,x0,v0,opt);
% figure(23);plot(x10(end-1,:),x10(end,:));

%% Visualization of the Results
disp('Visualization of the Results');

%% Plotting the profile of a single Periodic orbit
% If you have followed the GUI-tutorial, then uncomment the load lines, and
% comment out (x=xlc1/2);
% load Systems/JR_slow/diagram/LC_Part2.mat
x =xlc2;
ii=10; % The i-th point on the branch
jj=1;  % The j-th variable that we're going to plot.
ndim=8;% The number of state variables
ntst=40;ncol=4;  %From the continuation
mesh=sort([reshape(repmat(f(1:ntst,ii),1,ncol)+...
     repmat((0:ncol-1)/ncol,ntst,1).*repmat(diff(f(1:ntst+1,ii)),1,ncol),ntst*ncol,1); 1]);
figure(26);plot(mesh*x(end-1,ii),x(jj:ndim:end-2,ii))
xlabel('Time [s]','Fontsize',14);ylabel(strcat('Variable ',num2str(jj)),'Fontsize',14);

%% Plotting Max/min ghost bifurcations
figure(27);clf(27);hold on;
% load Systems/JR_slow/diagram/LC_Part1.mat
x =xlc2;f=flc2;
out=x(3:8:end-2,:)-(x(5:8:end-2,:)+x(7:8:end-2,:))/2;
extrema=nan(6,size(out,2));
for i=1:size(out,2)
    if any(abs(f(42:49,i))<1)
        a = diff(out([1:end 1],i));
        PeakIndices = find(diff(sign(a)))+1;
        temp =out(PeakIndices,i);
        extrema(1:size(temp,1),i)=temp;
    end
end
% AMP=[max(extrema); min(extrema)];
plot(x(end,:),sort(extrema),'linestyle','none','Marker','.');
% load Systems/JR_slow/diagram/LC_Part2.mat
x=xlc3;f=flc3;
out=x(3:8:end-2,:)-(x(5:8:end-2,:)+x(7:8:end-2,:))/2;
extrema=nan(6,size(out,2));
for i=1:size(out,2)
    if any(abs(f(42:49,i))<1)
        a = diff(out([1:end 1],i));
        PeakIndices = find(diff(sign(a)))+1;
        temp =out(PeakIndices,i);
        extrema(1:size(temp,1),i)=temp;
    end
end
plot(x(end,:),sort(extrema),'linestyle','none','Marker','.');
xlim([100 200]);
xlabel('II','Fontsize',14);ylabel('Extrema in signal','Fontsize',14);