% TIME-DELAY EMBEDDING ON El Nino/ASZIMMETRIKUS LORENZ
% BALADTA, 2012

% CREDIT GOES TO BÉLA WEISS FOR THE TIME-DEALY EMBEDDING FUNCTIONS
clear all

%%%%%%%% PARAMÉTEREK %%%%%%%%
sysdim = 2*length(el_nino_de);
% y0 = 0.1*ones(sysdim,1);
y0 = [-0.6 0  0.3 ];
B=102;
C=3; 
p=0.83;
params = [B C p];
dT = 0.01; % lépésköz
refresh_T = 0.5; % ennyit oldjunk meg egyszerre: refresh_T/dT
hist_len_f = 100; % ennyi látszódjon a múltból: refresh_T/dt*hist_len_f adat
N = 400; % hány itráción keresztül menjen a szimuláció: N*refresh_T/dt adatpont
% Y = nan(refresh_T/dT*hist_len_f+1,sysdim); 
% Y(:,1) = y0;
% t = nan(refresh_T/dT*hist_len_f+1,1);

correlation_threshold = 1/2.71; % 1/exp(1)
m = 3; % hány dimenzióban végezzük a TDE-t
sdim = 1; % melyik változó szerint végezzük a TDE-t

% diff. egyenket megoldó beállításai
options = odeset('RelTol',1e-3,'AbsTol',1e-6); % defaults
% options = odeset('RelTol',1e-6,'AbsTol',1e-7);

%%%%%%%% SZÁMOLÁS %%%%%%%%

[t,Y] = ode45(@(t,y) el_nino_de(t,y,params),0:dT:(N+1)*refresh_T,y0,options);
if false 
    figure
    hold on
    plot([-t(end:-1:2); t], xcorr(Y(:,sdim), 'coeff'))
    line([-t(end); t(end)],[correlation_threshold; correlation_threshold], 'LineWidth', 2, 'Color', 'r')
    xlabel('Time (s)'); ylabel('Autocorrelation')
end

%%
lag = correlation_function_thresholding(xcorr(Y(:,sdim), 'coeff'), correlation_threshold, [-t(end:-1:2); t]);
% figure; plot(Y./repmat(max(Y),length(Y),1)) % alapján a lag nem ugyanaz mindegyik dimenzióra, vagyis az semi mindegy, melyik dimenzió alapján végezzük a TDE-t
lag = 15;
YTDE = time_delay_embedding(Y(:,sdim),m,lag);



%%%%%%%% ÁBRÁZOLÁS %%%%%%%%
doimgwrite = false; % írjuk ki az animáció? (lassít)
figure;
set(gcf,'Color','w')

% ORIGINAL
subplot(7,2,1)
hcs11 = plot(t(1:refresh_T/dT),Y(1:refresh_T/dT,1));
ylabel('x(t)')
subplot(7,2,3)
hcs12 = plot(t(1:refresh_T/dT),Y(1:refresh_T/dT,2));
ylabel('y(t)')
subplot(7,2,5)
hcs13 = plot(t(1:refresh_T/dT),Y(1:refresh_T/dT,3));
xlabel('Time'); ylabel('z(t)')

% TDE
subplot(7,2,2)
hcs21 = plot(t(1:refresh_T/dT),YTDE(1:refresh_T/dT,1));
ylabel('x_{TDE}(t)')
subplot(7,2,4)
hcs22 = plot(t(1:refresh_T/dT),YTDE(1:refresh_T/dT,2));
ylabel('y_{TDE}(t)')
subplot(7,2,6)
hcs23 = plot(t(1:refresh_T/dT),Y(1:refresh_T/dT,3));
xlabel('Time'); ylabel('z_{TDE}(t)')

% Fázis
% subplot(7,2,7:10)
% hph = plot(Y(:,1),Y(:,2),'b',Y(:,4),Y(:,5),'r');
subplot(7,2,[7 9 11 13])
% hph(1) = plot(Y(1:refresh_T/dT,1),Y(1:refresh_T/dT,2),'b');
hph(1) = plot3(Y(1:refresh_T/dT,1),Y(1:refresh_T/dT,2),Y(1:refresh_T/dT,3),'b'); grid on
xlabel('x(t)'); ylabel('y(t)'); zlabel('z(t)')
view(16,18)
axis([-20 20 -1 1 -0.6 0.8])
subplot(7,2,[8 10 12 14])
% hph(2) = plot(YTDE(1:refresh_T/dT,1),Y(1:refresh_T/dT,2),'b');
hph(2) = plot3(YTDE(1:refresh_T/dT,1),YTDE(1:refresh_T/dT,2),YTDE(1:refresh_T/dT,3),'b'); grid on
xlabel('x_{TDE}(t)'); ylabel('y_{TDE}(t)'); zlabel('z_{TDE}(t)')
axis([-20 20 -20 20 -15 15])
view(16,18)

% Szinkronizáció
% subplot(7,2,11:14)
% hsyn = plot(Y(:,1),Y(:,4),'b',Y(:,2),Y(:,5),'r');
% subplot(7,2,[11 13])
% hsyn(1) = plot(Y(1:refresh_T/dT,2),YTDE(1:refresh_T/dT,2),'b'); 
% xlabel('y(t)'); ylabel('y_{TDE}(t)')
% subplot(7,2,[12 14])
% hsyn(2) = plot(Y(1:refresh_T/dT,3),YTDE(1:refresh_T/dT,3),'b'); 
% xlabel('z(t)'); ylabel('z_{TDE}(t)')

xyz = 'xyz';
title_str = sprintf(...
    'Time-Delay Embedding of the Asymmetric Lorenz System (El Nino system with parameters of Tung)\nReconstructed from %s(t) in %d-dim with a time lag of %d ms',...
    xyz(sdim),m,t(lag+1)*1000);
hs = annotation('textbox',[0.1 0.94 0.8 0.055],'String',title_str,'HorizontalAlignment','center','LineStyle','none');
pos = get(gcf, 'Position'); 
sc = 1.5; % 1.5
set(gcf,'Position',[pos(1)*0.25 pos(2)*0.25 pos(3)*sc pos(4)*sc])

if doimgwrite
    first = 0;
    frame = getframe(gcf); 
    [mov, cmap] = rgb2ind(frame.cdata, 256, 'nodither'); 
    % mov = rgb2ind(frame.cdata, cmap, 'nodither'); 
    for k = 1:first
        mov(:,:,1,1+k) = mov(:,:,1,1);
    end
    imwrite(mov, cmap, 'el_nino_animation.gif', 'DelayTime', 0.05, 'LoopCount', inf);
else
    % waitforbuttonpress
    pause(0.01)
end


tic
index_beg = 1;
index_end = 1;
for n = 1:N
    if index_end <= hist_len_f*refresh_T/dT
        % index_beg = index_beg; % marad
        index_end = index_end+refresh_T/dT;
        set(get(hcs11,'Parent'),'XLim',[t(index_beg) t(index_end)])
        set(get(hcs12,'Parent'),'XLim',[t(index_beg) t(index_end)])
        set(get(hcs13,'Parent'),'XLim',[t(index_beg) t(index_end)])
        set(get(hcs21,'Parent'),'XLim',[t(index_beg) t(index_end)])
        set(get(hcs22,'Parent'),'XLim',[t(index_beg) t(index_end)])
        set(get(hcs23,'Parent'),'XLim',[t(index_beg) t(index_end)])
    else
        index_beg = index_beg+refresh_T/dT;
        index_end = index_end+refresh_T/dT;
        set(get(hcs11,'Parent'),'XLim',[t(index_beg) t(index_end)])
        set(get(hcs12,'Parent'),'XLim',[t(index_beg) t(index_end)])
        set(get(hcs13,'Parent'),'XLim',[t(index_beg) t(index_end)])
        set(get(hcs21,'Parent'),'XLim',[t(index_beg) t(index_end)])
        set(get(hcs22,'Parent'),'XLim',[t(index_beg) t(index_end)])
        set(get(hcs23,'Parent'),'XLim',[t(index_beg) t(index_end)])
    end
    if length(YTDE)<index_end % a TDE miatt YTDE rövidebb, mint Y
        break
    end
    
    % ORIGINAL
    set(hcs11,'XData',t(index_beg:index_end),'YData',Y(index_beg:index_end,1));
    set(hcs12,'XData',t(index_beg:index_end),'YData',Y(index_beg:index_end,2));
    set(hcs13,'XData',t(index_beg:index_end),'YData',Y(index_beg:index_end,3));
    % TDE
    set(hcs21,'XData',t(index_beg:index_end),'YData',YTDE(index_beg:index_end,1));
    set(hcs22,'XData',t(index_beg:index_end),'YData',YTDE(index_beg:index_end,2));
    set(hcs23,'XData',t(index_beg:index_end),'YData',YTDE(index_beg:index_end,3));
    % Fázis
    % set(hph(1),'XData',Y(index_beg:index_end,1),'YData',Y(index_beg:index_end,2));
    set(hph(1),'XData',Y(index_beg:index_end,1),'YData',Y(index_beg:index_end,2),'ZData',Y(index_beg:index_end,3));
    % set(hph(2),'XData',YTDE(index_beg:index_end,1),'YData',YTDE(index_beg:index_end,2));
    set(hph(2),'XData',YTDE(index_beg:index_end,1),'YData',YTDE(index_beg:index_end,2),'ZData',YTDE(index_beg:index_end,3));
    % Szinkronizáció
    % set(hsyn(1),'XData',Y(index_beg:index_end,2),'YData',YTDE(index_beg:index_end,2))
    % set(hsyn(2),'XData',Y(index_beg:index_end,3),'YData',YTDE(index_beg:index_end,3))
    
    if doimgwrite
        frame = getframe(gcf); 
        mov = rgb2ind(frame.cdata, cmap, 'nodither');
        imwrite(mov, cmap, 'el_nino_animation.gif', 'WriteMode', 'append', 'DelayTime', 0.05);
    else
        drawnow
        % waitforbuttonpress
        % pause(0.01)
    end
end
toc

if doimgwrite
    last = 5;
    for k = 1:last
        imwrite(mov, cmap, 'el_nino_animation.gif','WriteMode', 'append', 'DelayTime', 0.05);
    end
end
