clear all

a1=0.08;
a2=0.23;
a3=.185; 
m1=15;
m2=7; % for
m3=2.5;
eps=1; 
del=1;
m2=5;

a1=0.08;a2=0.23;a3=.2;m1=10;m2=4;m3=3.5; eps = 1;del=1; % Hal & Waltman values
s=.8;m1=m1*s;m2=m2*s;m3=m3*s;
%eps=.01;

parameters=[a1 a2 a3 m1 m2 m3 eps del];

atp=0.23;  %%%%% input upper limit of parameter a2
abt=0.08; %%%%% input  lower limit of parameter a2

PIterate=25; %%%%% input:  Maximal number of bisection iteration to find the parameter a2, which aligns the shooting global unstable manifold 
             %%%%%        to the glabal stable manifold
             
error=1e-6; %%%%% input: Stopage error

NumberOfRun=50; %%%%% input: Number of initial points of the local unstable manifold

hx=5e-04; %%%%% .000055 input: Minimal offset of unstable manifold points away from the saddle focus 
%hx=2e-9;

SBisection=25; %%%%% input: Number of bisection to align the shooting unstable manifold orbit to the z=equ(3) plane.

options=odeset('e','on','BDF', 'on','RelTol',1e-16,'AbsTol',[1e-16;1e-16;1e-16]);
options1=odeset('e','off','BDF', 'on','RelTol',1e-16,'AbsTol',[1e-16;1e-16;1e-16]);


tic
[LocalUnstablemanifold, LocalStablemanifold, GlobalUnstablemanifold, GlobalStablemanifold, a2par, slnkidx, shilnikovorbit, searchonz, srcherror, sorbiterror, shooterror]=searchshilnikovorbit(atp, abt, parameters, PIterate, error, NumberOfRun, hx, SBisection, options, options1);
toc

%save SWShilnikovOrbit10.mat atp abt PIterate error NumberOfRun hx SBisection options options1 parameters LocalUnstablemanifold LocalStablemanifold GlobalUnstablemanifold GlobalStablemanifold a2par slnkidx shilnikovorbit searchonz srcherror sorbiterror shooterror
%load SWShilnikovOrbit20.mat atp abt PIterate error NumberOfRun hx SBisection options options1 parameters LocalUnstablemanifold LocalStablemanifold GlobalUnstablemanifold GlobalStablemanifold a2par slnkidx shilnikovorbit searchonz srcherror sorbiterror shooterror
%parameters(end)=.5;
%return
%%
aa=abs(srcherror);
%aa=abs(sorbiterror);
bb=aa';
cc=sum(bb);
[~,kk]=min(cc);
[~,kk]=min(sorbiterror(:,2));

%kk=length(a2par);
Unstablemanifold=LocalUnstablemanifold(:,:,kk);
Stablemanifold=GlobalStablemanifold(:,:,kk);
parameters(2)=a2par(kk);

%load ShilnikovOrbit1.mat atp abt PIterate NumberOfRun NumOfBiSection hx NumOfPartition NumberOfIterate SBisection options parameters Unstablemanifold Stablemanifold slnk shilnikovinitial srcherror

a1=parameters(1);
a2=parameters(2);
a3=parameters(3);
m1=parameters(4);
m2=parameters(5);
m3=parameters(6);
eps=parameters(7);
del=parameters(8);

figure(6)
hold off

nullclines
grid on
hidden off
axis([0 1 0 1 0 1])
view(60,18)
axis off

%%
%plotunstablemanifold1

plotunstablemanifold2

%%
tend=800;
init=Stablemanifold;
    [~,solX]=ode15s('chemostatmodel',[0,100],init,options,[parameters equ]);
                strt=1;
                xx=solX(strt:end,1);
                yy=solX(strt:end,2);
                zz=solX(strt:end,3);
                
[~,idx]=min(sqrt((xx-equ(1)).^2+(yy-equ(2)).^2+(zz-equ(3)).^2)); 

plot3(xx(1:idx),zz(1:idx),yy(1:idx),'r-','linewidth',3,'markersize',5)
plot3(init(1),init(3),init(2),'ro','linewidth',2,'markersize',10,'MarkerFaceColor',[1 1 1])
%%
pzz=psectionfamily(:,3)-Stablemanifold(3);
pyy=psectionfamily(:,2)-Stablemanifold(2);
pxx=psectionfamily(:,1)-Stablemanifold(1);
%dist=sqrt(pxx.^2+pyy.^2+pzz.^2);
%[~,slnk]=min(dist);
[aa,idx]=sort(pzz);
bb=pyy(idx);
cc=pxx(idx);
d=0;

plot3(cc(1:end-d)+Stablemanifold(1),aa(1:end-d)+Stablemanifold(3),bb(1:end-d)+Stablemanifold(2),'m-','linewidth',2) 


%%
% Plot the Shilnikov orbit
%

[t,solX]=ode15s('chemostatmodelus',[0,300],shilnikovorbit(kk,2:end),options1,parameters);
  
[~,tend]=min(abs(t-shilnikovorbit(kk,1)));

plot3(solX(1:tend,1),solX(1:tend,3),solX(1:tend,2),'k-','linewidth',2,'markersize',5)

sorbiterror(kk,:)

%%
k=length(a2par);
%sht=zeros(k,1);
figure(9)
hold off
d=0;
clr=colormap(jet(k));
for i=1:k%length(a2par)
pxx=GlobalUnstablemanifold(1:end-d,1,i)-GlobalStablemanifold(1,1,i);
pzz=GlobalUnstablemanifold(1:end-d,3,i)-GlobalStablemanifold(1,3,i);
pyy=GlobalUnstablemanifold(1:end-d,2,i)-GlobalStablemanifold(1,2,i);    
[aa,idx]=sort(pzz);
%[aa,idx]=sort(pxx);
bb=pyy(idx);
%plot(aa(1:end-d),bb(1:end-d),'g-',[0],[0],'ro','markersize',5,'linewidth',1) 
plot(aa(1:end-d),bb(1:end-d),'color',clr(i,:),'linewidth',2)    
hold on
%sht(i)=interp1(aa,bb,0,'linear');
end

plot([0],[0],'ro','markersize',10,'linewidth',2,'MarkerFaceColor',[1 1 1])    

%return
%%

figure(8)
hold off

nullclines
plot3(solX(:,1),solX(:,3),solX(:,2),'r-','linewidth',1,'markersize',5)

grid on
hidden off
set(gca,'fontsize',12)

xlabel('X','fontsize',16)
ylabel('Z','fontsize',16)
zlabel('Y','fontsize',16)

axis([0 1 0 1 0 1])
view(60,18)

return
%%

figure(6)

text(-.05,-.05,-.05, '0','fontsize',24)
text(.7,-.05,-.05, 'X','fontsize',24)
text(0,1,0.05, 'Z','fontsize',24)
text(-.05,-.05,.8, 'Y','fontsize',24)

axis off

% zoom-in axis([.69 .71 .158 .17 .104 .12])
%%
figure(7)

axis([-.02 0.02 -.02 .02])

ax = gca;
set(ax,'YTick', [-0.02 -0.01 0 0.01 0.02]);

set(gca,'fontsize',16)

xlabel('Z','fontsize',20)
ylabel('Y','fontsize',20)

%%
plot3([0,.5],[1,0],[0,0],'k--')
plot3([0,.5],[0,0],[1,0],'k--')