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

clear all
close all

%Run Time
tf = 400*(2*pi); %integer multiple of 2pi gives the number intersections with the poincare section.
%tf=17;
%tf =15;

%parameters
%d =0.25; g = 0.3;%d = 0.25
%d = 0.25; g = 0.1 % periodic orbit,

d = 0.15; g = 0.3; %%strange attactor


 %d = 0.22; g = 0.3; %%period 3 (implies chaos)
%


% d = 0.25; g = 0.1;%stable periodic orbit
 %


w = 1.0;
p = [d;g;w];

% %Equilibria
% xeq = a/5;
% yeq = 1+xeq^2;
% Ueq = [xeq;yeq];



%Initial condition(s)
U1 =[1,0];
%U2 =[0,3]; 
%U1 = X1 + 0.01*V_u;%[-0.4165 -0.9091] %0];%[0 1 0]
%U2 = X1-0.01*V_u;%[0 -1 0];  %%Trajectories on both sides of the unstable manifold Ws(C0)


% Generate Solution


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

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

% Create 2d Plot of solutions
figure(1)
plot(y1(:,1),y1(:,2))%,'linewidth',3)
%hold on
%lot(y2(:,1),y2(:,2),'linewidth',3)
%plot(xeq,yeq,'.k','MarkerSize',10)
%hold off
xlabel('x')
ylabel('y')

%%Create 3d plot, 
figure(2)
plot3(t1,y1(:,1),y1(:,2))%,'LineWidth',2)


figure(3) %omit transient)
plot(y1(50000:end,1),y1(50000:end,2))%,'linewidth',3)
%hold on
%lot(y2(:,1),y2(:,2),'linewidth',3)
%plot(xeq,yeq,'.k','MarkerSize',10)
%hold off
xlabel('x')
ylabel('y')

%%Construct Poincare map
dt = t1(2) - t1(1);
 Ip = find(abs(mod(t1,2*pi))<dt);%finds times when t is close to integer multiples of 2*pi
 xp = y1(Ip,1);
 yp = y1(Ip,2);
 figure(4)
 plot(xp,yp,'.','MarkerSize',10)
 
 
 figure(5) %omit transient
 xp1 = y1(Ip(125:end),1);
 yp1 = y1(Ip(125:end),2);
 plot(xp1,yp1,'.','MarkerSize',10)
 





