  function P2CVL2()
% Orbit Counting  around both attractors

% Parameters...
sigma = 10; rho = 28; beta = 8/3; eta = sqrt(beta*(rho-1));
u = [rho-1;eta;eta];           % Attractor 1
v = [rho-1;-eta;-eta];         % Attractor 2

% Our idea is to work with two zones. The comet at y is in Zone_u if
%                       (y-u)'*(v-u) <= 0
% and in Zone_v if
%                       (y-v)'*(u-v) <= 0
%
% If the comet is in Zone_u then it is ``behind'' Attractor 1 as seen 
% from an observer at Attractor 2.
% If the comet is in Zone_v then it is ``behind'' Attractor 2 as seen 
% from an observer at Attractor 1.
% Note that you cannot be in Zone_u and Zone_v at the same time.
%
% We define four events:
% 
%       Event 1 : Enter Zone_u
%       Event 2 : Leave Zone_u
%       Event 3 : Enter Zone_v
%       Event 4 : Leave Zone_v
% 
% Possible transitions and what we define them to mean:
%
%       Event 2 --> Event 1     Start of another orbit around Attractor 1
%       Event 2 --> Event 3     Move from Attractor 1 to Attractor 2.
%       Event 4 --> Event 3     Start of another orbit around Attractor 2
%       Event 4 --> Event 1     Move from Attractor 2 to Attractor 1.

% Simulate from t=0 to t = tFinal where
tFinal = 100;
% Initial position of the comet..
y0 = (u+v)/2 + 5*(1-.5*rand(3,1));



opts = odeset('RelTol',.000001,'events',@LorenzEvents);
[tOut,yOut,tE,yE,iE]= ode45(@Lorenz,[0 tFinal],y0,opts);
clc
% Tabulate the orbit history

orbVec = [];
OrbitCount=1;
m = length(iE)-1;
for k=1:m-1
    % Look at the transition from event iE(k) to event  iE(k+1)...
    if iE(k)==2 && iE(k+1)==1
        % Going around Attractor 1 again.
        OrbitCount = OrbitCount+1;
        
    elseif iE(k)==2 && iE(k+1)==3
        % Leaving Attractor 1 for Attractor 2
        % Record the number of orbits during this Attractor 1 phase..
        orbVec = [orbVec OrbitCount];
        % Initialize counter in anticipation of the upcoming Attractor 2
        % phase.
        OrbitCount = 1;
        
    elseif iE(k)==4 && iE(k+1)==3
        % Going around Attractor 2 again...
        OrbitCount = OrbitCount+1;
        
    elseif iE(k)==4 && iE(k+1)==1
        % Leaving Attractor 2 for Attractor 1
        % Record the number of orbits during this Attractor 2 phase..
        orbVec = [orbVec OrbitCount];
        % Initialize counter in anticipation of the upcoming Attractor 2
        % phase.
        OrbitCount = 1;
    end
end
fprintf('\n\nThe last transition: Event %1d --> Event %1d',iE(m-1),iE(m))
fprintf('\n\nTerminating value of OrbitCount = %1d\n',OrbitCount)
% If you want to watch the simulation...
close all
%shg
%comet3(yOut(:,1),yOut(:,2),yOut(:,3))
x = 1:length(orbVec);
h = .7;
bar(x(1),orbVec(1),h,'r')
hold on
for k=2:length(orbVec)
    if rem(k,2)==0
        bar(x(k),orbVec(k),h,'b');
    else
        bar(x(k),orbVec(k),h,'r');
    end
end
hold off
axis([0 length(orbVec)+1 0 max(orbVec)+1])
if iE(1)==1
   title('A1 Orbits in Red      A2 orbits in Blue','fontsize',16) 
else
   title('A1 Orbits in Blue      A2 orbits in Red','fontsize',16)
end
Total_A1_Orbits = sum(orbVec(1:2:length(orbVec)));
Total_A2_Orbits = sum(orbVec(2:2:length(orbVec)));
xlabel(sprintf('Total A1 Orbits = %3d   Total A2 Orbits = %3d',Total_A1_Orbits,Total_A2_Orbits),'fontsize',14)




  
function yp = Lorenz(t,y)
sigma = 10; rho = 28; beta = 8/3; 
yp = [-beta*y(1)+y(2)*y(3)    ;...
      -sigma*y(2)+sigma*y(3)  ; ...
      -y(2)*y(1)+rho*y(2)-y(3)...
      ];
  
  
 function [gstop,isterminal,direction] = LorenzEvents(t,y)
% Event function for the Lorenz simulation with
sigma = 10; rho = 28; beta = 8/3; eta = sqrt(beta*(rho-1));
u = [rho-1;eta;eta];           % Attractor 1
v = [rho-1;-eta;-eta];         % Attractor 2

% Comet is in Zone_u if c1<=0 where
c1 = (y-u)'*(v-u);
% Comet is in Zone_v if c2<=0 where
c2 = (y-v)'*(u-v);
% Event definitions:
%    1  Comet entering Zone_u
%    2  Comet leaving  Zone_u
%    3  Comet entering Zone_v
%    4  Comet leaving  Zone_v
gstop      = [ c1 c1 c2 c2];
direction  = [ -1 1  -1  1];
isterminal = [ 0 0 0 0];
