% ---------------------------------------------%
% LORENTZ DEMO PROGRAM 2016
% BY SVEND BERTEL DAHL-JENSEN
% SOLVES AND VISULAIZES THE SOLUTIONS OF THE LORENTZ EQUATIONS IN THE
% XZ PLANE.
% MAKES TWO INTIAL POINTS, ONE WITH A SMALL OFFSET, AND SHOWS THE SOLUTION
% TRJECTORIES LIVE TO SHOW THE EFFECT OF STRANGE ATTRACTORS AND CHAOS.
% HOPE YOU ENJOY THE CODE!
%-----------------------------------------------%

clear all
close all
clc


% HANDLES TO CHANGE 
live = 1; % 1 = LIVE IMAGING , 0 = ONLY END PLOT
tail = 10; % LENGTH OF SOLUTION TAIL DURING LIVE
offset = [0 eps 0]; % OFFSET OF THE TWO INITAL POINTS, [] FOR 1 SOLUTION 
% LORENTZ PARAMETERS, SEE CHAPTER 9.5 STROGATZ 2006 FOR A OVERVIEW OF THE
% PARAMETER SPACE
sigma=10;
beta=8/3;
rho=28;
% TIME PARAMETERS
t_min = 0;
t_max = 100;
% INTIAL POINTS STARTING COORDINATES
% EXTEND INIT ARRAY FOR MULTIPLE TRAJECTORIES (BEYOND THE WHAT IS GAINED WITH OFFSET)
x_init= [0];
y_init= [1];
z_init= [0];

%'f' IS THE SET OF DIFFERENTIAL EQUATIONS AND 'a' IS A VECTOR CONTAINING THE THREE PARAMETERS
f = @(t,a) [-sigma*a(1) + sigma*a(2); rho*a(1) - a(2) - a(1)*a(3); -beta*a(3) + a(1)*a(2)];


figure(1)
set(1,'Color',[1 1 1])
axis([-30 30 0 60])
for list = 1:length(x_init)
    % SOLVES 'f' WITH RUNGE KUTTA
    [ts,ys] = ode45(f,[t_min,t_max],[x_init(list),y_init(list),z_init(list)]);
    if isempty(offset) == 0
    [ts2,ys2] = ode45(f,[t_min,t_max],[x_init(list)+offset(1),y_init(list)+offset(2),z_init(list)+offset(3)]);
    end
    % MAKES LIVE ANIMATION
    if live ==1
        if isempty(offset) == 0
            for t_int=t_min+1:size(ys,1)
                plot(gca,ys((max(t_int-tail,1):t_int),1),ys((max(t_int-tail,1):t_int),3),'Tag','live_line')
                hold on
                plot(gca,ys2((max(t_int-tail,1):t_int),1),ys2((max(t_int-tail,1):t_int),3),'Tag','live_line','Color',[1 0 0])
                hold off
                axis([-30 30 0 60])
                %axis off
                drawnow
                delete(findall(gca,'Tag','live_line'))
            end
        else
            for t_int=t_min+1:size(ys,1)
                plot(gca,ys((max(t_int-tail,1):t_int),1),ys((max(t_int-tail,1):t_int),3),'Tag','live_line')
                axis([-30 30 0 60])
                drawnow
                delete(findall(gca,'Tag','live_line'))
            end
        end
    end
    % PLOTS END RESULT
    hold on
    plot(ys(:,1),ys(:,3))
    plot(ys2(:,1),ys2(:,3),'r-')
    plot(ys(1,1),ys(1,3),'*','Color',[0 0 0]) % START POINT
    plot(ys(end,1),ys(end,3),'g*') % ENDING POINT
    hold off    
end
hold off
legend('Vector field','Trajectory','Start point','End point')



