%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         9. practice         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nonlinear dynamics systems  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         2016.11.16.         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    Ljapunov, bifurcation    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

kezdet=0.2; %0.3; %arbitrary initial conditions;
iteraciok=2000;
mu=0.1:0.001:4;% the full diagram
%mu=3.8:0.00001:3.9; %periodic window 

% Ricker iteration
%a=0.1;
%mu=0.1:0.01:30;

Population=zeros(length(mu), iteraciok); %store the data here: rows: the mu parameters, columns: the iterations
Population(:,1)=kezdet*ones(length(mu),1); %initial condition for every mu

for i=1:(iteraciok-1) %do the iteration
    Population(:, i+1)=mu'.*Population(:, i).*(1-Population(:, i));
    %Population(:,i+1)=mu'.*Population(:,i).*exp(-a*Population(:,i)); %Ricker

end

%%
size(repmat(mu',[1 101]))   % only show the last 100 element
size(Population(:,900:1000))

Population_derivalt = repmat(mu',[1 101])-2.*repmat(mu',[1 101]).*Population(:,900:1000);
%Population_derivalt = repmat(mu',[1 101]).*exp(-a*Population(:,900:1000)).*(-a*Population(:,900:1000)+1); %Ricker

Ljapunov = (log(abs(prod(Population_derivalt,2))))/100;
%Ljapunov = sum(log(abs(Population_derivalt)),2)/100;

%%

figure;
subplot(2,1,1);
ax(1)=subplot(2,1,1);
%plot(mu, Population(:, 0.99*iteraciok:iteraciok), 'b.', 'MarkerSize', 1);
plot(mu, Population(:, 0.99*iteraciok:iteraciok), 'b.', 'MarkerSize', 1,'Parent',ax(1));
title('Bifurkációs diagram');
xlabel('mu');
ylabel('populáció');


subplot(2,1,2);
ax(2)=subplot(2,1,2);
plot(mu,Ljapunov);
hold on;
%plot([min(mu),max(mu)],[0,0],'r')
plot([min(mu),max(mu)],[0,0],'r','Parent',ax(2));
title('Ljapunov');
xlabel('mu');
ylabel('Ljapunov exponens');

linkaxes(ax,'x');
