% Illustrate sensitivity to initial conditions.

global global_beta = 8/3
global global_sigma = 10
global global_rho = 28

times = 0:0.01:50; 

% Fixed points of the Lorenz system
Origin = [0,0,0]
Cplus = [sqrt(global_beta*(global_rho-1)), \
         sqrt(global_beta*(global_rho-1)), \
         global_rho-1]
Cminus = [-sqrt(global_beta*(global_rho-1)), \
         -sqrt(global_beta*(global_rho-1)), \
         global_rho-1]

% Initial conditions. Chosen to be near the attractor but visible on plot.
initial =  [-5,10,30]

trajectory1 = lsode("lorenz", initial, times);

% Second trajectory with randomly perturbed initial condition.
% The perturbation, epsilon, is a vector of fixed length and random
% direction.
epsilon = randn(1,3); % normal_rnd(0,1,1,3);
epsilon = 1e-10*epsilon/norm(epsilon);
trajectory2 = lsode("lorenz", initial+epsilon, times);

% Analyze the difference between the trajectories

displacement = trajectory2 - trajectory1;
distance = sqrt(sum(displacement.^2, 2));


%Create distance plot.
figure(1)
hold off
semilogy(times, distance)
title("Sensitivity to Initial Conditions")
xlabel("t")
ylabel("| trajectory1 - trajectory2 |")
print("distance.pdf",'-dpdfwrite')


% Plots
figure(2)
hold off  % Let's start from scratch.

elevation = 23;
azimuth = 20;
x1 = trajectory1(:,1);
y1 = trajectory1(:,2);
z1 = trajectory1(:,3);
plot3(x1,y1,z1,'b','LineWidth',5)
   % Seems like Octave doesn't pick up the linewidth...
hold on
x2 = trajectory2(:,1);
y2 = trajectory2(:,2);
z2 = trajectory2(:,3);
plot3(x2,y2,z2,'r','LineWidth',5)
title('Lorenz Attractor: Two Trajectories')

plot3(initial(1), initial(2), initial(3), 'k*');
plot3(Origin(1),  Origin(2),  Origin(3),  'go');
plot3(Cplus(1),   Cplus(2),   Cplus(3),   'g+');
plot3(Cminus(1),  Cminus(2),  Cminus(3),  'gx');
legend('Trajectory 1', 'Trajectory 2',...
   'Initial condition', 'Origin', 'C+', 'C-')
view(azimuth, elevation)
xlabel('x')
ylabel('y')
zlabel('z')
print('attractor.pdf','-dpdfwrite')

%Output frames for a movie of the attractor.
for frame = 1:360;
      azimuth = mod(frame + 30, 360);
         % view doesn't work with azimuth > 360. Lame!
      view(azimuth, elevation);
      filename = strcat('/tmp/movie', sprintf('%04u',frame), '.png')
      print(filename,'-dpng','-S1024,768');
end

