%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% solves the Lorenz equations with ode45        %
% requires the use of file:  myode.m            %
% produces two figures:                         %
% figure 1 = x,y,z vs. time                     %
% figure 2 = 3D trajectory of "lorenz butterfly %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function lorenz_kktung

clear;
% declare global parameters 
global r;
global sigma;
global b;

% set parameters 
r = 28;
sigma = 10;
b = 8/3;

% initial conditions
x0 = 3;
y0 = 15;
z0 = 1;
v0 = [x0,y0,z0];
t0 = 0;
tf = 60;
tspan = [t0 tf];

% solves 'myode' using ode45 routine

[t,vout] = ode45(@myode,tspan,v0);

% note:
paramstring = sprintf( ', {\\itr}=%d, {\\it\\sigma}=%d,{\\itb}=%2.2f', r, sigma, b);
icstring = sprintf( ', ({\\itx}_0,{\\ity}_0,{\\itz}_0=(%d,%d,%d)', x0, y0, z0 );

figure(1)
clf

%plots x(t) in the top 1/3 of the figure
subplot(3,1,1)
plot( t, vout(:,1) );
xlabel( '{\itt}' );
ylabel( '{\itx}({\itt})' );
title( ['Lorenz Solutions', paramstring, icstring ] )
axis tight;
grid on;

%plots y (t) in the second 1/3 of the figure
subplot(3,1,2)
plot( t, vout(:,2) );
xlabel( '{\itt}' );
ylabel( '{\ity}({\itt})' );
axis tight;
grid on;

%plots z (t) in the third 1/3 of the figure
subplot(3,1,3)
plot( t, vout(:,3) );
xlabel( '{\itt}' );
ylabel( '{\itz}({\itt})' );
axis tight;
grid on;

figure(2);
clf;
% plots trajectory 
plot3( vout(:,1), vout(:,2), vout(:,3) );
hold on;
% plot a circle at initial point
plot3( vout(1,1), vout(1,2), vout(1,3), 'o' );

axis tight;
xlabel( '{\itx}' );
ylabel( '{\ity}' );
zlabel( '{\itz}' );
title( ['Lorenz Trajectory', paramstring, icstring] );

camproj('perspective');
campos([180.8444 -309.0346 241.0763] );
camtarget( [1.6213   3.2901   25.8520] );

end

%-------------------

function vprime = myode( t, v )
% parameters (defined elsewhere)
global sigma;
global r;
global b;

% map current vector to user-friendly variable names
x = v(1);
y = v(2);
z = v(3);

% compute derivatives with the Lorenz model
xprime = sigma*(-x+y);
yprime = r*x-y-x*z;
zprime = -b*z + x*y;

% map back to vector format
vprime=zeros(3,1);
vprime(1) = xprime;
vprime(2) = yprime;
vprime(3) = zprime;

end