%% Numerical Continuation for Unstable Hopf-Cycle.   
clc
clear all
%Parameter values for plotting switching manifold 
alpha=0.624; er=1; vreset=0.1538; sjump=0.8; ts=1.4; wjump = 0.0189;
tw = 65; vpeak=1.4615; er =1;
%Bifurcation parameters 
I = 0.4; 
g = 1.2308

%Run an initial simulation to determine the equilibrium point
int = [0,0]
[t,y] = ode45(@(t,y) IzDIRECT(t,y,I,g),[0,200],int');
init
xeq = y(end,:)
[x0,v0]=init_EP_EP(@Izmeanfield2,xeq',[I;g],[1])

%Configure Numerical Constants for continuation
opt=contset;
  % opt=contset(opt,'Adapt',10000);
opt=contset(opt,'Backward',1);
opt=contset(opt,'MaxStepSize',0.001);
opt=contset(opt,'MinStepSize',0.00001); 
   %opt=contset(opt,'MaxNumPoints',10000);
opt=contset(opt,'Singularities',2);
% opt=contset(opt,'Increment',1e-5); 
opt =contset(opt,'Eigenvalues',1);

%Equilibrium continuation 
  [x,v,s,h,f]=cont(@equilibrium,x0,[],opt);
  
%Equilibrium point, Parameter at Hopf Bifurcation Point 
x1 = x(1:2,s(2).index); p =[x(end,s(2).index);g];
[x0,v0]=init_H_LC(@Izmeanfield2,x1,p,[1],1e-6,20,4);
     
%Change some numerical constants/parameters for numerical continuation
%of a limit cylce
opt = contset(opt,'MaxNumPoints',50);
opt = contset(opt,'Multipliers',1);
opt=contset(opt,'MaxStepSize',100);
opt=contset(opt,'MinStepSize',0.1); 
opt=contset(opt,'Adapt',50);

%Numerical continuation of a limit cycle, and plotting results.  Plots the
%unstable limit cylce in Blue
[xlc,vlc,slc,hlc,flc]=cont(@limitcycle,x0,v0,opt);
xlc(end,:) = xlc(end,:)*2.5*65*65;  x(end,:) = x(end,:)*2.5*65*65; %(dimensionalization, I_app = I * k*vr*vr
plotcycle(xlc,vlc,slc,[size(xlc,1) 1 2]), hold on
cpl(x,v,s,[size(x,1),1,2])
    
 %% Direct Simulations for Bursting Limit Cycle
 I1min = min(xlc(end,:))/(2.5*65*65); %Range of I values
 I1max = max(xlc(end,:))/(2.5*65*65);
 dI = (I1max-I1min)/(20);  %step over the range

 I1 = I1min; %Initialize at min and increase in increments of DI 
tspan = 0:1:100; %Times to compute the limit cycle.  
index = 0; 
REC1 = []; REC2 = []; REC3 = []; 
  while I1<I1max 
      index = index + 1; 
 I1 = I1 + dI ;
 [t,y] = ode45(@(t,y) IzDIRECT(t,y,I1,g),[0,200],zeros(2,1)); %Run direct simulations, get rid of initial transient
 ynot = y(end,:);
 [t,y] = ode45(@(t,y) IzDIRECT(t,y,I1,g),tspan,ynot'); %Stable Limit Cycle 
 %drawnow
 if mod(index,2) == 1 
 plot3(I1*2.5*65*65 + 0*y(:,1),y(:,1),y(:,2),'k'), hold on %Plot the limit cycle for fixed I in 3 dimensions.  
 end
 %We need to save these to generate the bursting limit cycle surface.  
REC1(:,index) =  I1*2.5*65*65 + 0*y(:,1);
REC2(:,index) = y(:,1);
REC3(:,index) = y(:,2); 

  end
  %Plot the switching manifold in RED 
  axis([I1min*2.5*65*65,I1max*2.5*65*65+20,0,0.3,0.05,0.15])
  s = 0.2*(0:0.1:1); u = 0.2*(0:0.1:1);
  s = ones(11,1)*s; u = ones(11,1)*u;
  I = -2.5*65*65*(-u' + g*er*s - 0.25*(alpha+g*s).^2);
  surf(I,s,u','FaceColor','red','EdgeColor','none')
  clear alpha 
colormap hsv


%Plot the stable bursting limit cycle manifold in Green 
surf(REC1,REC2,REC3,'FaceColor','green','EdgeColor','none')
  clear alpha 
colormap hsv
alpha(.2)
% for i = 1:10
%     plot3(REC1(10*i,:),REC2(10*i,:),REC3(10*i,:),'k'), hold on 
% end
xlabel('$I_{app}$','Interpreter','LateX','FontSize',14)
 ylabel('$s$','Interpreter','LateX','FontSize',14)
 zlabel('$\langle w \rangle$','Interpreter','LateX','FontSize',14)