
% Code to extract a Poincare section
% RCG, Fa'20

%Initial conditions/setup
ttrans=100;
torbit=100;
tpoincare=1000;
IC0=[8,3,33];
params=[10,28,8/3]; % params [sigma,r,b]
options1 = odeset('Reltol',1.e-5);
options2 = odeset('Reltol',1.e-5,'Refine',3,'Events',@myevent);
ms=10; % MarkerSize
plotfrac=1; % fraction of orbit to plot

%Just integrate to get rid of transients:
[Ttemp,Ytemp]=ode45(@(t,y) LorenzEqs(t,y,params),[0,ttrans],IC0,options1); 
IC=Ytemp(end,:); % this eliminates transients

%Now integrate to save orbit:
[To,Yo]=ode45(@(t,y) LorenzEqs(t,y,params),[0,torbit],IC,options1); 

%Now the Poincare section
[T,Y,Te,Ye,Ie]=ode45(@(t,y) LorenzEqs(t,y,params),[0,tpoincare],IC,options2);

%Plot IC, orbit, and PS
iin=[1:round(length(To)*plotfrac)]; % plot just a fraction of data
figure(1);clf
plot3(Yo(iin,1),Yo(iin,2),Yo(iin,3))
hold on
plot3(Yo(1,1),Yo(1,2),Yo(1,3),'r.','MarkerSize',ms*2)
plot3(Ye(:,1),Ye(:,2),Ye(:,3),'k.','MarkerSize',ms)
hold off
xlabel('x'); ylabel('y'); zlabel('z');
view([22 27])

%Plot just the PS
figure(2);clf
plot(Ye(:,2),Ye(:,3),'k.','MarkerSize',ms)
xlabel('y'); ylabel('z');title('Poincare section')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [value,isterminal,direction] = myevent(t,y) % event function
  CUT = 0;
  value = y(1)-CUT; % Poincare section when y(1) = CUT
  isterminal = 0; % don't quit when first zero is found
  direction = 1; % -1,+1 or 0 for both
end

