% main

% Initial conditions
t0 = 0;
w0 = [1,1,1];

% Maximum time of integration
T = 100;

% Step definitions
h0 = 0.01;              % initial step size
h = h0;                 % used step size
hmin = (1e-5)*h0;       % minimum step size allowed

% Model parameters and definitions
sigma = 10;
rho = 28;
beta = 8/3;

F = @(s,w) Lorenz(s,w,sigma,rho,beta);

% Poincare Section definition
Section = @ (w) w(3) - 15;
tsec = [];
wsec = [];

% Result allocation
t = zeros(100000,1);
w = zeros(100000,max(size(w0)));

t(1) = t0;
w(1,:) = w0;
counter  = 2;

% Initiate figure object
figure;
az = 0; el = 0;     % viewpoint definition

%%%%%% ITERATE
while t(counter) <= T
    
    % iterate a step
    [t(counter), w(counter,:)] = RK4(t(counter-1),w(counter-1,:),h,F); 
    
    % evaluete if solutuion passes through the poincare section
    if or(sign(Section(w(counter,:))) > sign(Section(w(counter-1,:))), sign(Section(w(counter,:))) < sign(Section(w(counter-1,:))))
        h = h/2;
        if h <= hmin
           h = h0;
           tsec = [tsec;t(counter)];
           wsec = [wsec;w(counter,:)];
           counter = counter + 1;
        end
    else
        
        % plot result
        plot3(w(1:counter,1),w(1:counter,2),w(1:counter,3))
        xlabel('x'); ylabel('y'); zlabel('z')
        hold on;
        % plot points that pass through Poincare section
        if max(size(wsec))>0
            plot3(wsec(:,1),wsec(:,2),wsec(:,3),'o');
        end
        % plot current point in trajectory
        plot3(w(counter,1),w(counter,2),w(counter,3),'*')
        hold off
        % rotate and view solution
        az = mod(az + 0.5,360);
        view(az,el)
        pause(0.005);
        
        %iterate
        counter = counter + 1;
    end
    
    
end

% eliminate unused entries in solution vector 
t(counter+1:end) = [];
w(counter+1,:) = [];