clear

%% parameter values

% system parameters
epsilon=0.05;
gamma=0.97;
mu2=0.12;
alpha3=0.3;
beta2=0;
beta3=0;

num_sampl=1000;
mu1=linspace(0.02,0.07,num_sampl);

% range of initial conditions for x0 and x0d
x10_range=[1, 2];
x1d0_range=[1, 2];

% time span and time interval
tfinal=1000;
t_steps=1000;
tspan=linspace(0,tfinal,t_steps);

%% generate random initial conditions

x0=zeros(num_sampl,4);
x0(:,1)=rand(num_sampl,1)*(x10_range(2)-x10_range(1))+x10_range(1);
x0(:,2)=rand(num_sampl,1)*(x1d0_range(2)-x1d0_range(1))+x1d0_range(1);
x0(:,3)=x0(:,1);


%% perform simulations
options = odeset('RelTol', 1e-9, 'AbsTol', 1e-9);
x=cell(length(mu1),1);
for i=1:num_sampl
    [~, x{i}] = ode45(@(t,x) system_VDP(t, x, epsilon, gamma, mu1(i), mu2, alpha3), tspan, x0(i,:), options);
end

%% save data
clear("i");
save VDP_4incremental_testing_raw_simulation_data.mat

%% normalize data polar_01_log_mvm

clear

load VDP_4incremental_testing_raw_simulation_data.mat

rho1=cell(length(x),1);
rho2=cell(length(x),1);
for i=1:length(x)
    polar_pure1=sqrt(x{i}(:,1).^2+x{i}(:,2).^2);
    polar_pure2=sqrt(x{i}(:,3).^2+x{i}(:,4).^2);
    max_pol1=max(polar_pure1);
    max_pol2=max(polar_pure2);
    polar_011=polar_pure1/max_pol1;
    polar_012=polar_pure2/max_pol2;
    polar_01_log1=log(polar_011);
    polar_01_log2=log(polar_012);
    for j=20:length(polar_01_log1)
        if max([polar_01_log1(j-19:j)])<-14
            polar_01_log1=polar_01_log1(1:j);
            break
        end
    end
    polar_01_log_resumpled1=spline(linspace(1,length(polar_01_log1),length(polar_01_log1)), polar_01_log1, linspace(1,length(polar_01_log1),1000));
    rho1{i}=movmean(polar_01_log_resumpled1,round(length(polar_01_log_resumpled1)/20));
    for j=20:length(polar_01_log2)
        if max([polar_01_log2(j-19:j)])<-14
            polar_01_log2=polar_01_log2(1:j);
            break
        end
    end
    polar_01_log_resumpled2=spline(linspace(1,length(polar_01_log2),length(polar_01_log2)), polar_01_log2, linspace(1,length(polar_01_log2),1000));
    rho2{i}=movmean(polar_01_log_resumpled2,round(length(polar_01_log_resumpled2)/20));
end


% save normalized data
clear("x");
clear("i");
save VDP_4incremental_testing_polar_01_log_mvm.mat



%% normalize data minmax [-1 1]

clear

load VDP_4incremental_testing_raw_simulation_data.mat

rho1=cell(length(x),1);
rho2=cell(length(x),1);
for i=1:length(x)
    massimo1=max(x{i}(:,1));
    minimo1=min(x{i}(:,1));
    minmax1=zeros(length(x{i}(:,1)),1);
    massimo2=max(x{i}(:,3));
    minimo2=min(x{i}(:,3));
    minmax2=zeros(length(x{i}(:,1)),1);
    for j=1:length(x{i})
        if x{i}(j,1)>0
            minmax1(j)=x{i}(j,1)/massimo1;
        else
            minmax1(j)=x{i}(j,1)/abs(minimo1);
        end
        if x{i}(j,1)>0
            minmax2(j)=x{i}(j,3)/massimo2;
        else
            minmax2(j)=x{i}(j,3)/abs(minimo2);
        end
    end
    for j=50:length(minmax1)
        if max([minmax1(j-49:j)])<1e-14
            minmax1=minmax1(1:j);
            break
        end
    end
    for j=50:length(minmax2)
        if max([minmax2(j-49:j)])<1e-14
            minmax2=minmax2(1:j);
            break
        end
    end
    minmax_resumpled1=spline(linspace(1,length(minmax1),length(minmax1)), minmax1, linspace(1,length(minmax1),1000));
    minmax_resumpled2=spline(linspace(1,length(minmax2),length(minmax2)), minmax2, linspace(1,length(minmax2),1000));
    rho1{i}=minmax_resumpled1;
    rho2{i}=minmax_resumpled2;
end


% save normalized data
clear("x");
clear("i");
save VDP_4incremental_testing_minmax.mat

%% normalize data polar_01

clear

load VDP_4incremental_testing_raw_simulation_data.mat

rho1=cell(length(x),1);
rho2=cell(length(x),1);
for i=1:length(x)
    polar_pure1=sqrt(x{i}(:,1).^2+x{i}(:,2).^2);
    polar_pure2=sqrt(x{i}(:,3).^2+x{i}(:,4).^2);
    max_pol1=max(polar_pure1);
    max_pol2=max(polar_pure2);
    polar_011=polar_pure1/max_pol1;
    polar_012=polar_pure2/max_pol2;
    for j=20:length(polar_011)
        if max([polar_011(j-19:j)])<-14
            polar_011=polar_011(1:j);
            break
        end
    end
    polar_01_resumpled1=spline(linspace(1,length(polar_011),length(polar_011)), polar_011, linspace(1,length(polar_011),1000));
    rho1{i}=polar_01_resumpled1;
    for j=20:length(polar_012)
        if max([polar_012(j-19:j)])<-14
            polar_012=polar_012(1:j);
            break
        end
    end
    polar_01_resumpled2=spline(linspace(1,length(polar_012),length(polar_012)), polar_012, linspace(1,length(polar_012),1000));
    rho2{i}=polar_01_resumpled2;
end


% save normalized data
clear("x");
clear("i");
save VDP_4incremental_testing_polar_01.mat





%% normalize data polar_01_mvm

clear

load VDP_4incremental_testing_raw_simulation_data.mat

rho1=cell(length(x),1);
rho2=cell(length(x),1);
for i=1:length(x)
    polar_pure1=sqrt(x{i}(:,1).^2+x{i}(:,2).^2);
    polar_pure2=sqrt(x{i}(:,3).^2+x{i}(:,4).^2);
    max_pol1=max(polar_pure1);
    max_pol2=max(polar_pure2);
    polar_011=polar_pure1/max_pol1;
    polar_012=polar_pure2/max_pol2;
    for j=20:length(polar_011)
        if max([polar_011(j-19:j)])<-14
            polar_011=polar_011(1:j);
            break
        end
    end
    polar_01_resumpled1=spline(linspace(1,length(polar_011),length(polar_011)), polar_011, linspace(1,length(polar_011),1000));
    rho1{i}=movmean(polar_01_resumpled1,round(length(polar_01_resumpled1)/20));
    for j=20:length(polar_012)
        if max([polar_012(j-19:j)])<-14
            polar_012=polar_012(1:j);
            break
        end
    end
    polar_01_resumpled2=spline(linspace(1,length(polar_012),length(polar_012)), polar_012, linspace(1,length(polar_012),1000));
    rho2{i}=movmean(polar_01_resumpled2,round(length(polar_01_resumpled2)/20));
end


% save normalized data
clear("x");
clear("i");
save VDP_4incremental_testing_polar_01_mvm.mat


%% normalize data polar_01_log

clear

load VDP_4incremental_testing_raw_simulation_data.mat

rho1=cell(length(x),1);
rho2=cell(length(x),1);
for i=1:length(x)
    polar_pure1=sqrt(x{i}(:,1).^2+x{i}(:,2).^2);
    polar_pure2=sqrt(x{i}(:,3).^2+x{i}(:,4).^2);
    max_pol1=max(polar_pure1);
    max_pol2=max(polar_pure2);
    polar_011=polar_pure1/max_pol1;
    polar_012=polar_pure2/max_pol2;
    polar_01_log1=log(polar_011);
    polar_01_log2=log(polar_012);
    for j=20:length(polar_01_log1)
        if max([polar_01_log1(j-19:j)])<-14
            polar_01_log1=polar_01_log1(1:j);
            break
        end
    end
    polar_01_log_resumpled1=spline(linspace(1,length(polar_01_log1),length(polar_01_log1)), polar_01_log1, linspace(1,length(polar_01_log1),1000));
    rho1{i}=polar_01_log_resumpled1;
    for j=20:length(polar_01_log2)
        if max([polar_01_log2(j-19:j)])<-14
            polar_01_log2=polar_01_log2(1:j);
            break
        end
    end
    polar_01_log_resumpled2=spline(linspace(1,length(polar_01_log2),length(polar_01_log2)), polar_01_log2, linspace(1,length(polar_01_log2),1000));
    rho2{i}=polar_01_log_resumpled2;
end


% save normalized data
clear("x");
clear("i");
save VDP_4incremental_testing_polar_01_log.mat





















