%%    UMD AOSC 615   SPRING 2009
%%    Lorenz 3 variable model ('63)

%% startup
clear all
close all
cpu0=cputime;
disp(sprintf('\n lorenz63: %s',datestr(now)))

%% global
global ode_prmtr

%% input
% system parameters
a=10; 
b=8/3; 
r=28;
% inital condition (before spin-up)
x0=[8; 8; 3];
% integration time
Tsu=10;   %  initial spin-up
Tot=20;       %  integration time
dTot=0.1;   %  output
% plot
c_txyz=['b','r','g'];
axs_txyz=[0 Tot -30 50];
axs_xyz=[-30 30 -30 30 0 60];

%% setup
% dim
N=length(x0);            
% output 
tot=[0:dTot:Tot];
KTot=length(tot);
xot=zeros(N,KTot);
% model parameter for ode
ode_prmtr=[a b r];    

%% initial spinup
[tsu,xsu]=ode45(@lorenz63_dxdt,[0 Tsu],x0);
xot(:,1)=xsu(end,:)'; 

%% model forward integration
for k=2:KTot;  kn1=k-1; xk0=xot(:,kn1);
    [tk,xk]=ode45(@lorenz63_dxdt,[tot(kn1) tot(k)],xk0);
    xot(:,k)=xk(end,:)'; 
end

%% plot
% t and xk
figure
for n=1:N
    plot(tot,xot(n,:),c_txyz(n),'LineWidth',2); if(ishold~=1); hold on; end
end
legend('x_1','x_2','x_3','Location','SouthEast')
axis(axs_txyz)
xlabel('time'); ylabel('x_1, x_2,  x_3')

% 3d
figure
plot3(xot(1,:),xot(2,:),xot(3,:),'b','LineWidth',1)
xlabel('x_1'),ylabel('x_2'),zlabel('x_3')
axis(axs_xyz);

%% finish
disp(sprintf(' ... end of lorenz63 (cpu=%s sec)\n',num2str(cputime-cpu0)))
