%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Solve the Lorenz 3-D System.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

           function lorenz

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  We define the program as a function
%  to enable the computation of the
%  right-hand side to be called as
%  a subfunction (in the same file).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all    % Clear all variables from memory

global r     % Parameters to be passed to sub-function
global sigma %
global b     %

r = 28;      % Define parameter values for system
sigma = 10;  % 
b = 8/3;     %

x0 = 1; y0 = 1; z0 = 1;

xyz0 = [x0; y0; z0];   %  Array of initial values

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Flush out the first 200 values
%  to eliminate any transients.
%
t0 = 0; 
tfinal = 200;
tspan = [t0 tfinal];
[time xyz] = ode45(@ROC_Lorenz,tspan,xyz0);
nt = size(time,1);
x0 = xyz(nt,1); y0 = xyz(nt,2); z0 = xyz(nt,3);

xyz0 = [x0; y0; z0];   %   Re-set initial values

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Solve over a span for the "steady" solution.
t0 = 0; 
tfinal = 100;
tspan = [t0 tfinal];
[time xyz] = ode45(@ROC_Lorenz,tspan,xyz0);

nt = size(time,1);
ns = 1;
x = xyz (ns:nt,1);
y = xyz (ns:nt,2);
z = xyz (ns:nt,3);
t = time(ns:nt);

%   Plot the solution.
subplot(3,2,1); plot(t,x)
title('(a) x versus t')
subplot(3,2,3); plot(t,y)
title('(b) y versus t')
subplot(3,2,5); plot(t,z)
title('(c) z versus t')

subplot(3,2,2); plot(x,y); axis('square')
title('(d) x versus y')
subplot(3,2,4); plot(y,z); axis('square')
title('(e) y versus z')
subplot(3,2,6); plot(x,z); axis('square')
title('(f) x versus z')

%  Save the image
print -deps 'Lorenz1.eps'
pause

clf
plot(x,z); axis('square'); grid on
title('\bf THE LORENZ ATTRACTOR')
C1x = +sqrt(b*(r-1)); C1y = C1x; C1z = r-1;
C2x = -sqrt(b*(r-1)); C2y = C2x; C2z = r-1;
text(C1x,C1z,'\bf .C1')
text(C2x,C2z,'\bf .C2')
print -deps 'Lorenz2.eps'
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf
plot3(x,y,z); %  axis('square'); grid on
xlabel('x'); ylabel('y'); zlabel('z');
title('\bf THE LORENZ ATTRACTOR')
%%% print -deps 'Lorenz2.eps'
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Plot twin solutions.

nt = size(time,1);
x01 = xyz(nt,1); y01 = xyz(nt,2); z01 = xyz(nt,3);
x02 = x01*(1+0.004);y02 = y01; z02 = z01;

xyz01 = [x01; y01; z01];
xyz02 = [x02; y02; z02];

t0 = 0; 
tfinal = 12;
tspan = [t0 tfinal];
[t1 xyz1] = ode45(@ROC_Lorenz,tspan,xyz01);
[t2 xyz2] = ode45(@ROC_Lorenz,tspan,xyz02);

x1 = xyz1(:,1); y1 = xyz1(:,2); z1 = xyz1(:,3);
x2 = xyz2(:,1); y2 = xyz2(:,2); z2 = xyz2(:,3);

clf
plot(t1,x1,'b','LineWidth',2); grid on
hold on
plot(t2,x2,'--r','LineWidth',2);
hold off
title('\bf IDENTICAL TWIN EXPERIMENT: x(t)')

%  Save the image
print -djpeg 'Lorenz3.jpg'
print -deps 'Lorenz3.eps'

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   rhs of Lorenz 3-D system
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xyzdot = ROC_Lorenz(t,xyz)

global r
global sigma
global b

x = xyz(1); y = xyz(2); z = xyz(3);

xdot = sigma*(y-x);
ydot = r*x - y - x*z;
zdot = x*y-b*z;

xyzdot = [xdot; ydot; zdot];

% size(xyzdot)
% pause

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
