%%%%Script to solve lorenz equations,
%%%% Change constants in function file lorenz.m
clear all
close all

%Run Time
tf = 300;
%tf=17;
%tf =15;


%parameters
sigma = 10;
b = 8/3;
%r = 24;%heteroclinic orbits 2-13
%r = 13.9; %close to homoclinic
%r=13.926;%homoclinic value
rhb = (sigma*(sigma+b+3))/(sigma-b-1);
%r = rhb-0.4;%%hopf value
r = 28;  %Lorenz's fateful parameters


%%Other cool stuff
%r=21; %transient chaos
%r = 24.5; %chaos and stable fixed point


%r = 280;%200; 200-300 interesting periodic orbits
%r = 147.5;%period doubling values 350,313,260,224,222,218,216.2, 147.5,
%r = 216.2
%r = 218
%r = 222
%r = 224
%%r = 260
%r = 313
%r = 350
%%for these don't plot the initial transient t\in [0 10]


p = [sigma,r,b];

%%Equilibria
C0 = [0,0,0]';
Cm = [-sqrt(b*(r-1)),-sqrt(b*(r-1)),r-1]';
Cp = [+sqrt(b*(r-1)),+sqrt(b*(r-1)),r-1]';


%%Jacobian Matrix
xs = C0(1);
ys = C0(2);
J = [-sigma sigma 0 ; r -1 -xs;ys xs -b];
[VV,DD] = eig(J)
[ds,Is] = sort(diag(DD));
%vus = VV(:,Is(3));
vus = VV(:,Is(3));
%%Jacobian Matrix
xs = Cm(1);
ys = Cm(2);
J = [-sigma sigma 0 ; r -1 -xs;ys xs -b];
[VV,DD] = eig(J)

%Initial condition(s)
U1 =0.0001*vus';%[-0.4165 -0.9091 0];%[0 1 0]
U2 = -U1;%[0 -1 0];  %%Trajectories on both sides of the unstable manifold Ws(C0)

%U1 = [10 .1 60];%[40 -4 20]; %for transient chaos
%U2 = -U1;

%%bistability with chaos
%U1 = Cp+[1,1,1]';
%U2 = Cm+[1,1,1]';

% Generate Solution


ff =@(tt,yy) lorenz(tt,yy,p);

opts = odeset('MaxStep',2e-3);
%%Solutions for different initial conditions
[t1, y1] = ode15s(ff, [0 tf], U1,opts);
[t2, y2] = ode15s(ff, [0 tf], U2,opts);

% Create 3d Plot
figure(1)
plot3(y1(:,1), y1(:,2), y1(:,3));
hold on
plot3(y2(:,1), y2(:,2), y2(:,3));
plot3(C0(1),C0(2),C0(3),'.k','MarkerSize',10)
plot3(Cp(1),Cp(2),Cp(3),'.k','MarkerSize',10)
plot3(Cm(1),Cm(2),Cm(3),'.k','MarkerSize',10)
hold off
xlabel('x')
ylabel('y')
zlabel('z')


% Create 3d Plot - without transient
IT = floor(length(t1)/2) : length(t1);
figure(11)
plot3(y1(IT,1), y1(IT,2), y1(IT,3));
hold on
plot3(y2(IT,1), y2(IT,2), y2(IT,3));
plot3(C0(1),C0(2),C0(3),'.k','MarkerSize',10)
plot3(Cp(1),Cp(2),Cp(3),'.k','MarkerSize',10)
plot3(Cm(1),Cm(2),Cm(3),'.k','MarkerSize',10)
hold off
xlabel('x')
ylabel('y')
zlabel('z')

% % Create 3d Plot
% figure(1)
% plot3(y1(5000:end,1), y1(5000:end,2), y1(5000:end,3));
% hold on
% plot3(y2(5000:end,1), y2(5000:end,2), y2(5000:end,3));
% plot3(C0(1),C0(2),C0(3),'.k','MarkerSize',10)
% plot3(Cp(1),Cp(2),Cp(3),'.k','MarkerSize',10)
% plot3(Cm(1),Cm(2),Cm(3),'.k','MarkerSize',10)
% hold off
% xlabel('x')
% ylabel('y')
% zlabel('z')




%%Plot times series in x and z
figure(2)
subplot(2,1,1)
plot(t1,y1(:,1))
xlabel('t')
ylabel('x')
subplot(2,1,2)
plot(t1,y1(:,3))
xlabel('t')
ylabel('z')
