clear
clear functions
% Solve and Plot the Lorenz System of 3 Differential Equations
% initial data
x0=[10 10 10];

% initialize initial time and end time
t0=0; tf=100; 
[xout, tout]=bs_lorenz(x0, t0, tf);

% repeat Lorenz plot with slith perturbation in the initial conditions
x0(1)=x0(1)+0.2;
[xout2, tout2]=bs_lorenz(x0, t0, tf);



% analyze the phase space
% plot the solution curve and set various display options
figure(1); 
clf  % clear figure content
hold on
% Plot output of first initial condition set
hp=plot3(xout(:,1), xout(:,2), xout(:,3),'color','b'); 
plot3(xout(end,1), xout(end,2), xout(end,3),'*g'); 
% Plot output of second initial condition set
hp=plot3(xout2(:,1), xout2(:,2), xout2(:,3),'color','r'); 
plot3(xout2(1,1), xout2(1,2), xout2(1,3),'*r'); 
% change some of the graph properties
set(hp,'LineWidth',0.5);   
box on;   
xlabel('x','FontSize',14);
ylabel('y','FontSize',14);   
zlabel('z','FontSize',14);
title('Lorenz Attractor phase space','FontSize',14);
axis([-20 20 -40 40 0 60]);
set(gca,'CameraPosition',[200 -200 200],'FontSize',14);
set(hp,'LineWidth',0.5);



% analyze the time series of each variable
figure(2)
clf
variab='XYZ';
for ivar=1:3
subplot(3,1,ivar)
hold on
plot(tout,xout(:,ivar),'b'); 
plot(tout2,xout2(:,ivar),'r');
title ( variab(ivar) )
end


% plot the PDFs of X, Y and Z
figure(3)
clf
variab='XYZ';
for ivar=1:3
subplot(3,1,ivar)
hold on
[n, edges]=bs_sample_pdf_1d(xout(:,ivar));
bar(edges,n,'histc');
title ( ['PDF of ',variab(ivar)] );
end


% 2D view of phase space in X Z
figure(4); 
subplot(2,1,1)
hold on
hp=plot(xout(:,1), xout(:,3),'color','b'); 
plot(xout(1,1), xout(1,3),'*b'); 
% change some of the graph properties
set(hp,'LineWidth',0.5);   
box on;   
xlabel('X','FontSize',14);
ylabel('Z','FontSize',14);   
title('Phase diagram for X and Z','FontSize',14);   
axis([-20 20 0 60]);
% repeat Lorenz plot with slith perturbation in the initial conditions
hp=plot(xout2(:,1), xout2(:,3),'color','r'); 
plot(xout2(1,1), xout2(1,3),'*r'); 
set(hp,'LineWidth',0.5);
subplot(2,1,2)
% Plot Joint PDF of X Z
tmp1=[xout(:,1); xout2(:,1)];
tmp2=[xout(:,3); xout2(:,3)];
[jpdf, axis1, axis2]=bs_sample_pdf_2d(tmp1,tmp2);
nlev=60;
utl_contourfill(axis1,axis2,jpdf,nlev); 
m=colormap; m(1:4,:)=1; colormap(m); colorbar
title('Joint PDF of X and Z','FontSize',14);   
xlabel('X','FontSize',14);
ylabel('Z','FontSize',14);   

% 2D view of phase space in X Z for Y close to 0
figure(5); 
subplot(2,1,1)
range=[-1 1];
hold on
in=find( xout(:,2) < range(end) & xout(:,2) > range(1) );
hp=plot(xout(in,1), xout(in,3),'color','b'); 
% change some of the graph properties
set(hp,'LineWidth',0.5);   
box on;   
xlabel('X','FontSize',14);
ylabel('Z','FontSize',14);   
title('Phase diagram for X and Z (when -1 < Y < 1)','FontSize',14);   
axis([-20 20 0 60]);
tmp1=[xout(in,1)];
tmp2=[xout(in,3)];
in=find(xout2(:,2) < range(end) & xout2(:,2) > range(1) );
hp=plot(xout2(in,1), xout2(in,3),'color','r'); 
set(hp,'LineWidth',0.5);
tmp1=[tmp1; xout2(in,2)];
tmp2=[tmp2; xout2(in,3)];
subplot(2,1,2)
% Plot Joint PDF of X Z for given Y interval
[jpdf, axis1, axis2]=bs_sample_pdf_2d(tmp1,tmp2);
nlev=60;
utl_contourfill(axis1,axis2,jpdf,nlev); 
m=colormap; m(1:4,:)=1; colormap(m); colorbar
title('Conditional Joint PDF of X and Z (for -1 < Y < 1)','FontSize',14);   
xlabel('X','FontSize',14);
ylabel('Z','FontSize',14);