%simualte cnn
%close all;
clear all;
clc;

options = odeset('RelTol',1e-5, 'AbsTol', 1e-10, 'MaxStep', 1e-1);

%Max_time = 100;
% for lady's shoe:
Max_time = 1500;

A=[ 0  1.2 0;...
    0  2   0;...
    0 -1.2 0];
B=[ 0 0 0;...
    0 1 0;...
    0 0  0];
P.Z=0;

P.A=reshape(A,1,9);
P.B=reshape(B,1,9);

P.Boundary='Constant';
P.BounValue=0;
P.SizeX=2;
P.SizeY=1;
P.Input = @(t)([4.04*sin(pi*t/2); 0]);
InitState= [0.1 0.1]; %[0.1 0.2];

% define every point of tspan (and not just the first and last one), this way 
% the solutions will be provided over the specified domain-point
freq = 100;
[t,X] = ode45(@cell_equation_boundary,[0:(1/freq):Max_time-1/freq],InitState,options,P);

%out=reshape(X,size(t,1),P.SizeX,P.SizeY);

figure;
hold on;
plot(t, X(:,1), 'r.-');
plot(t, X(:,2), 'b.-');
xlabel('time (t)', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('state (x)', 'FontSize', 12, 'FontWeight', 'bold');
title('x1 and x2 in time domain; chaotic behaviour', 'FontSize', 14);

figure;
hold on;
plot(X(:, 1), X(:, 2), 'b.-');
xlabel('x1', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('x2', 'FontSize', 12, 'FontWeight', 'bold');
title('phase-space portrait; chaotic behaviour', 'FontSize', 14);

%%{
% for the next figure, better to have
% Max_time at least 1500
lady_index = (mod(t, 4)>=3.97 | mod(t, 4)<=0.03); %(mod(t, 4)==0);
t_lady = t(lady_index);
figure;
plot(X(lady_index, 1), X(lady_index, 2), 'b.');
xlabel('x1', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('x2', 'FontSize', 12, 'FontWeight', 'bold');
title('sampled phase-space; chaotic behaviour', 'FontSize', 14);
%%}




% compute the frequency spectra
m = size(X, 1);                 % Window length
n = pow2(nextpow2(m));          % Transform length
y1 = fft(squeeze(X(:, 1)), n);  % DFT
y2 = fft(squeeze(X(:, 2)), n);  % DFT
f = (0:n-1)*(freq/n);           % Frequency range
power1 = y1.*conj(y1)/n;        % Power of the DFT
power2 = y2.*conj(y2)/n;        % Power of the DFT

figure;
subplot(1, 2, 1);
hold on;
plot(f(1:floor(n/2)), power1(1:floor(n/2)), 'r.-');
plot(f(1:floor(n/2)), power2(1:floor(n/2)), 'b.-');
xlabel('Frequency (Hz)', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('Power', 'FontSize', 12, 'FontWeight', 'bold');
title('Spectrum', 'FontSize', 14);

subplot(1, 2, 2);
hold on;
plot(f(1:floor(n/2)), power1(1:floor(n/2)), 'r.-');
plot(f(1:floor(n/2)), power2(1:floor(n/2)), 'b.-');
xlim([0 1]);
xlabel('Frequency (Hz)', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('Power', 'FontSize', 12, 'FontWeight', 'bold');
title('Spectrum, enlarged', 'FontSize', 14);


