% Delay coordinate reconstruction of Lorenz 
% attractor from x(t) timeseries
% for apm147
% (dsa, 11/28/04)
function doall=recon_lorenz_phase_space
clear all;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% integrate system
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma=10.0;
b=8.0/3.0;
r=28;


% integration time:
T=60;
%set error bounds for integration
options=odeset('RelTol',10^(-6),'AbsTol',[10^(-6)*ones(1,3)]);
% integrate:
[t,x] = ode45(@lor,[0 T],[7.95 7.95 23.74],options,[sigma; b ;r]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% use delay coordinate reconstruction
% reconstruct phase space from x time series
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% interpolate onto even time grid
% (ode45 changes time-stepping)
dt=0.005;
tp=(0:dt:T)';
xp=interp1(t,x(:,1),tp);
NP=length(tp);

tau=0.15;
y=zeros((NP-2*tau/dt),3);
y(:,1) = xp(1:(NP-2*tau/dt));
y(:,2) = xp((tau/dt+1):(NP-tau/dt));
y(:,3) = xp((2*tau/dt+1):NP);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
subplot(2,1,1)
plot3(x(:,1),x(:,2),x(:,3))
xlabel('x')
ylabel('y')
zlabel('z')
title('Phase Space')
view(75,15)
subplot(2,1,2)
plot3(y(:,1),y(:,2),y(:,3))
xlabel('x(t)')
ylabel('x(t-\tau)')
zlabel('x(t-2\tau)')
title(['Reconstructed Phase Space, \tau = ' num2str(tau)])
view(75,15)

figure(2)
subplot(2,1,1)
plot(x(:,1),x(:,3))
xlabel('x')
ylabel('z')
title('Phase Space')
subplot(2,1,2)
plot(y(:,1),y(:,2))
xlabel('x(t)')
ylabel('x(t-\tau)')
title(['Reconstructed Phase Space, \tau = ' num2str(tau)])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy = lor(t,x,kh)
sigma=kh(1); b=kh(2); r=kh(3);
dy=[sigma*(x(2)-x(1)),r*x(1)-x(2)-x(1)*x(3),x(1)*x(2)-b*x(3)]';
