clear

%% parameter values

parameters_PnP;

ucr=0.91107;

num_sampl=1000;
u=linspace(ucr*0.9,ucr*1.1,num_sampl);

% range of initial conditions 
h0_range=[0.02, 0.1]/10;
hd0_range=[0.02, 0.05]/10;
a0_range=[0.2, 0.6]/2;
ad0_range=[0.2, 0.6]/2;

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

%% generate random initial conditions

x0(:,1)=rand(num_sampl,1)*(h0_range(2)-h0_range(1))+h0_range(1);
x0(:,3)=rand(num_sampl,1)*(hd0_range(2)-hd0_range(1))+hd0_range(1);
x0(:,2)=rand(num_sampl,1)*(a0_range(2)-a0_range(1))+a0_range(1);
x0(:,4)=rand(num_sampl,1)*(ad0_range(2)-ad0_range(1))+ad0_range(1);


%% perform simulations
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-8);
x=cell(length(u),1);
for i=1:num_sampl
    i
    par=coefficients_PnP(xa,ra,beta,nu,Omega,zetaa,zetah,xih3,xia3,xia5,u(i));
    [~, x{i}] = ode45 (@(t, xt) system_PnP(t,xt,par), tspan, x0(i,:), options);
end

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

%% normalize data polar_01_log_mvm

clear

load PnP_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}(:,3).^2);
    polar_pure2=sqrt(x{i}(:,2).^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)])<-7
            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)])<-7
            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 PnP_4incremental_testing_polar_01_log_mvm.mat


%% normalize data minmax [-1 1]

clear

load PnP_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}(:,2));
    minimo2=min(x{i}(:,2));
    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,2)/massimo2;
        else
            minmax2(j)=x{i}(j,2)/abs(minimo2);
        end
    end
    for j=50:length(minmax1)
        if max([minmax1(j-49:j)])<1e-7
            minmax1=minmax1(1:j);
            break
        end
    end
    for j=50:length(minmax2)
        if max([minmax2(j-49:j)])<1e-7
            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 PnP_4incremental_testing_minmax.mat


%% normalize data polar_01

clear

load PnP_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}(:,3).^2);
    polar_pure2=sqrt(x{i}(:,2).^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)])<-7
            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)])<-7
            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 PnP_4incremental_testing_polar_01.mat


%% normalize data polar_01_mvm

clear

load PnP_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}(:,3).^2);
    polar_pure2=sqrt(x{i}(:,2).^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)])<-7
            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)])<-7
            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 PnP_4incremental_testing_polar_01_mvm.mat


%% normalize data polar_01_log

clear

load PnP_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}(:,3).^2);
    polar_pure2=sqrt(x{i}(:,2).^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)])<-7
            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)])<-7
            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 PnP_4incremental_testing_polar_01_log.mat














