clear all

global cds sys

sys.gui.pausespecial=0;  %Pause at special points 
sys.gui.pausenever=1;    %Pause never 
sys.gui.pauseeachpoint=0; %Pause at each point

syshandle=@actinwavesys4again;  %Specify system file

SubFunHandles=feval(syshandle);  %Get function handles from system file
RHShandle=SubFunHandles{2};      %Get function handle for ODE


k0=.0; %Set ODE parameter
gam=1; %Set ODE parameter
k=.4; %Set ODE parameter - This is A0 in the actin waves paper
del=1; %Set ODE parameter
%hill=3; %This one is hardcoded
s1=.7 ;
s2=.7 ;
r0=.5 ;
kn=1 ;
ks=.25 ;
Tn=1; %Set ODE parameter
tau=.1 ;

xinit=[.1,.1,0,0]; %Set ODE initial condition

%Specify ODE function with ODE parameters set
RHS_no_param=@(t,x)RHShandle(t,x,k0,gam,k,del,s1,s2,r0,kn,ks,Tn,tau); 

%Set ODE integrator parameters.
options=odeset;
options=odeset(options,'RelTol',1e-5);
options=odeset(options,'maxstep',1e-1);

%Integrate until a steady state is found.
[tout xout]=ode45(RHS_no_param,[0,200],xinit,options);

%%
%%%%% Continuation from equilibrium %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Set initial condition as the endpoint of integration.  Use
%to bootstrap the continuation.
xinit=xout(size(xout,1),:);

pvec=[k0,gam,k,del,s1,s2,r0,kn,ks,Tn,tau]';      % Initialize parameter vector

ap=1;

[x0,v0]=init_EP_EP(syshandle, xinit', pvec, ap); %Initialize equilibrium


opt=contset;
opt=contset(opt,'MaxNumPoints',400); %Set numeber of continuation steps
opt=contset(opt,'MaxStepsize',.01);  %Set max step size
opt=contset(opt,'Singularities',1);  %Monitor singularities
opt=contset(opt,'Eigenvalues',1);    %Output eigenvalues 
opt = contset(opt,'InitStepsize',0.01); %Set Initial stepsize

[x1,v1,s1,h1,f1]=cont(@equilibrium,x0,v0,opt);
figure()
cpl(x1,v1,s1,[5 1]);

% %%
% %%%%% Continuation from equilibrium backward %%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [x0,v0]=init_EP_EP(syshandle, xinit', pvec, ap); %Initialize equilibrium
% opt=contset(opt,'Backward',1);
% [x2,v2,s2,h2,f2]=cont(@equilibrium,x0,v0,opt);
% % cpl(x2,v2,s2,[3 1]);

%%
%%%%% Branch swiching and continuation %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xBP=x1(1:4,s1(5).index);       %Extract branch point
pvec(ap)=x1(5,s1(5).index);  %Extract branch point 

%%%%%% Initialize branch point information %%%%%
[x0,vO]=init_BP_EP(syshandle, xBP, pvec, s1(5), 0.01);  

opt=contset(opt,'MaxNumPoints',150);
[x3,v3,s3,h3,f3]=cont(@equilibrium,x0,v0,opt); %Switch branches and continue.
figure()
cpl(x3,v3,s3,[5 1]);

%%
%%%%% Branch swiching and continuation %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% xBP=x1(1:2,s1(2).index);       %Extract branch point
% pvec(ap)=x1(3,s1(2).index);  %Extract branch point 
% 
% opt=contset(opt,'Backward',1);
% 
% %%%%%% Initialize branch point information %%%%%
% [x0,vO]=init_BP_EP(syshandle, xBP, pvec, s1(2), .01);  
% 
% [x4,v4,s4,h4,f4]=cont(@equilibrium,x0,v0,opt); %Switch branches and continue.


xNearBP=x3(1:4,10);
pvec(ap)=x3(5,10);

[x0,v0]=init_EP_EP(syshandle, xNearBP, pvec, ap); %Initialize equilibrium

opt=contset(opt,'MaxNumPoints',335);
opt=contset(opt,'Backward',1);
[x4,v4,s4,h4,f4]=cont(@equilibrium,x0,v0,opt); %Switch branches and continue.
figure()
cpl(x4,v4,s4,[5 1]);
%%
%%%%% Plotting script.  x=continuation info  f=eigenvalues %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



figure();
% colordef(gcf,'black')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xeqcurve=x1;
minevaleq=f1(4,:); %This is the last eigenvalue.  That is the one that determines stability

L=length(xeqcurve(1,:));

curveind=1;
lengthind=0;
maxlengthind=0;
evalstart=floor(heaviside(minevaleq(1)));
datamateq=zeros(4,L);

for i=1:L
    evalind=floor(heaviside(minevaleq(i)));
    if evalstart~=evalind
        curveind=curveind+1;
        i;
        evalstart=evalind;
        maxlengthind=max(lengthind,maxlengthind);   
        lengthind=0;
    end
    datamateq(1,i)=xeqcurve(5,i); % This is the parameter that is varied.
    datamateq(2,i)=xeqcurve(1,i); % This is the dependent axis of the bifurcation plot.  The one you wish to plot
    datamateq(3,i)=evalind;
    datamateq(4,i)=curveind;  
    
    lengthind=lengthind+1;
end

maxlengthind=max(maxlengthind,lengthind);

curveindeq=curveind;

for i=1:curveindeq
    index=find(datamateq(4,:)==i);
    eval(['curve' num2str(i) 'eq' '=datamateq(1:3,index);']);
end

for i=1:curveindeq
    stability=eval(['curve' num2str(i) 'eq(3,1)']);
    if stability==0
        plotsty='-';
    else
        plotsty=':';
    end
    
    plotcolor='k';

    plotstr=strcat(plotcolor,plotsty);
    
    plot(eval(['curve' num2str(i) 'eq(1,:)']),eval(['curve' num2str(i) 'eq(2,:)']),plotstr,'Linewidth',4)
    hold on
end


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% xeqcurve=x2;
% minevaleq=f2(2,:); %This is the last eigenvalue.  That is the one that determines stability
% 
% L=length(xeqcurve(1,:));
% 
% curveind=1;
% lengthind=0;
% maxlengthind=0;
% evalstart=floor(heaviside(minevaleq(1)));
% datamateq=zeros(4,L);
% 
% for i=1:L
%     evalind=floor(heaviside(minevaleq(i)));
%     if evalstart~=evalind
%         curveind=curveind+1;
%         i;
%         evalstart=evalind;
%         maxlengthind=max(lengthind,maxlengthind);   
%         lengthind=0;
%     end
%     datamateq(1,i)=xeqcurve(3,i); % This is the parameter that is varied.
%     datamateq(2,i)=xeqcurve(1,i); % This is the dependent axis of the bifurcation plot.  The one you wish to plot
%     datamateq(3,i)=evalind;
%     datamateq(4,i)=curveind;  
%     
%     lengthind=lengthind+1;
% end
% 
% maxlengthind=max(maxlengthind,lengthind);
% 
% curveindeq=curveind;
% 
% for i=1:curveindeq
%     index=find(datamateq(4,:)==i);
%     eval(['curve' num2str(i) 'eq' '=datamateq(1:3,index);']);
% end
% 
% for i=1:curveindeq
%     stability=eval(['curve' num2str(i) 'eq(3,1)']);
%     if stability==0
%         plotsty='-';
%     else
%         plotsty=':';
%     end
%     
%     plotcolor='k';
% 
%     plotstr=strcat(plotcolor,plotsty);
%     
%     plot(eval(['curve' num2str(i) 'eq(1,:)']),eval(['curve' num2str(i) 'eq(2,:)']),plotstr,'Linewidth',4)
%     hold on
% end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xeqcurve=x3;
minevaleq=f3(4,:); %This is the last eigenvalue.  That is the one that determines stability

L=length(xeqcurve(1,:));

curveind=1;
lengthind=0;
maxlengthind=0;
evalstart=floor(heaviside(minevaleq(1)));
datamateq=zeros(4,L);

for i=1:L
    evalind=floor(heaviside(minevaleq(i)));
    if evalstart~=evalind
        curveind=curveind+1;
        i;
        evalstart=evalind;
        maxlengthind=max(lengthind,maxlengthind);   
        lengthind=0;
    end
    datamateq(1,i)=xeqcurve(5,i); % This is the parameter that is varied.
    datamateq(2,i)=xeqcurve(1,i); % This is the dependent axis of the bifurcation plot.  The one you wish to plot
    datamateq(3,i)=evalind;
    datamateq(4,i)=curveind;  
    
    lengthind=lengthind+1;
end

maxlengthind=max(maxlengthind,lengthind);

curveindeq=curveind;

for i=1:curveindeq
    index=find(datamateq(4,:)==i);
    eval(['curve' num2str(i) 'eq' '=datamateq(1:3,index);']);
end

for i=1:curveindeq
    stability=eval(['curve' num2str(i) 'eq(3,1)']);
    if stability==0
        plotsty='-';
    else
        plotsty=':';
    end
    
    plotcolor='r';

    plotstr=strcat(plotcolor,plotsty);
    
    plot(eval(['curve' num2str(i) 'eq(1,:)']),eval(['curve' num2str(i) 'eq(2,:)']),plotstr,'Linewidth',2.5)
    hold on
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xeqcurve=x4;
minevaleq=f4(4,:); %This is the last eigenvalue.  That is the one that determines stability

L=length(xeqcurve(1,:));

curveind=1;
lengthind=0;
maxlengthind=0;
evalstart=floor(heaviside(minevaleq(1)));
datamateq=zeros(4,L);

for i=1:L
    evalind=floor(heaviside(minevaleq(i)));
    if evalstart~=evalind
        curveind=curveind+1;
        i;
        evalstart=evalind;
        maxlengthind=max(lengthind,maxlengthind);   
        lengthind=0;
    end
    datamateq(1,i)=xeqcurve(5,i); % This is the parameter that is varied.
    datamateq(2,i)=xeqcurve(1,i); % This is the dependent axis of the bifurcation plot.  The one you wish to plot
    datamateq(3,i)=evalind;
    datamateq(4,i)=curveind;  
    
    lengthind=lengthind+1;
end

maxlengthind=max(maxlengthind,lengthind);

curveindeq=curveind;

for i=1:curveindeq
    index=find(datamateq(4,:)==i);
    eval(['curve' num2str(i) 'eq' '=datamateq(1:3,index);']);
end

for i=1:curveindeq
    stability=eval(['curve' num2str(i) 'eq(3,1)']);
    if stability==0
        plotsty='-';
    else
        plotsty=':';
    end
    
    plotcolor='r';

    plotstr=strcat(plotcolor,plotsty);
    
    plot(eval(['curve' num2str(i) 'eq(1,:)']),eval(['curve' num2str(i) 'eq(2,:)']),plotstr,'Linewidth',2.5)
    hold on
end

axis([0 .5 0 .8])
scalefactor=3;
fontsizevar=11;
fontsizevar=fontsizevar*scalefactor;

xlabel('a','fontsize',fontsizevar,'color',[0 0 0])
ylabel('u^l','fontsize',fontsizevar,'color',[0 0 0])
title('Wave Pinning LPA','fontsize',fontsizevar,'color',[0 0 0])
set(gcf,'units','inches','pos',[0 0 3.1*scalefactor 2.61*scalefactor])
set(gca,'fontsize',fontsizevar,'YColor',[0 0 0],'XColor',[0 0 0],...
    'linewidth',1.5);

% set(gca,'fontsize',fontsizevar,'YColor',[.98 .98 .98],'XColor',[.98 .98 .98], ...
%     'Color',[0.87 0.92 0.98])
% axes('Parent',gca,'YColor',[1 1 1],'XColor',[1 1 1])


