clear all
close all
clc
set(0,'defaultTextInterpreter','latex');


%% FILE DESCRIPTION

% This file is used to assemble a population of sOCs starting from the
% ditribution of random sOC parameters.

% REFERENCES

% [1] L Romano, P Johannesson, F Bruzelius, and B Jacobson. " An enhanced
% stochastic operating cycle description including weather and traffic
% models". Transportation research part D: transport and environment 97,
% (2021).
% Available from: https://doi.org/10.1016/j.trd.2021.102878.

% [2] L Romano, P Johannesson, E Nordström, F Bruzelius, R Andersson, and B
% Jacobson. "A classification method of road transport missions and
% applications using the operating cycle format". IEEE Access 10, (2022).
% Available from: https://doi.org/10.1109/ACCESS.2022.3188872.

% [3] L Romano, M Godio, P Johannesson, F Bruzelius, T Ghandriz, and B
% Jacobson. "Development of an operating cycle for long-haul heavy-duty
% trucks: theory, analysis and applications". Submitted for publication.

% Please, cite the above references when using this script for purposes of
% academic or industrial research


%% INITIALISE PARAMETERS

% Number of sOCs to generate
dim = 10000; % Should be less or equal to the length of a variable in sOC

% Load the sOC variables
designationDir = 'C:\Users\luigir\Desktop\sOC for TRD\generated_sOCs';
sOC = load('C:\Users\luigir\Desktop\sOC for TRD\distributionsGenerated_SOC.mat');
sOC = sOC.sOC;
Ltot = sOC.Ltot*1000; % The mission length needs to be expressed in m here
Ls = 25;

% The parameter to the OC are the autoregressive coefficient for topography
phi_Y1 = sin(pi/2-2*Ls./sOC.L_h1);
phi_Y2 = sin(pi/2-2*Ls./sOC.L_h2);
phi_Y3 = sin(pi/2-2*Ls./sOC.L_h3);


%% ASSEMBLE POPULATION OF sOCs

% This should go from i=1:dim
for i = 1:1

    % Road types
    statOC.Road.TC.RoadTypes = {'City','Rural','Highway'};
    statOC.Road.TC.L = Ls;
    statOC.Road.TC.LL = [sOC.L_R1(i), sOC.L_R2(i), sOC.L_R3(i)];
    statOC.Road.TC.PP = [0 sOC.p_R12(i) sOC.p_R13(i);
        sOC.p_R21(i) 0 sOC.p_R23(i);
        sOC.p_R31(i) sOC.p_R32(i) 0];
    statOC.Road.TC.Ltot = Ltot(i);
    statOC.Road.TC.Speed = [50, 80, 100]/3.6;


    % Topography
    % Urban
    statOC.Road.Topo{1,1}.type = 'G1';
    statOC.Road.Topo{1,1}.a = phi_Y1(i);
    statOC.Road.Topo{1,1}.sy = sOC.Sigma_Y1(i);
    statOC.Road.Topo{1,1}.L = Ls;
    % Rural
    statOC.Road.Topo{1,2}.type = 'G1';
    statOC.Road.Topo{1,2}.a = phi_Y2(i);
    statOC.Road.Topo{1,2}.sy = sOC.Sigma_Y2(i);
    statOC.Road.Topo{1,2}.L = Ls;
    % Highway
    statOC.Road.Topo{1,3}.type = 'G1';
    statOC.Road.Topo{1,3}.a = phi_Y3(i);
    statOC.Road.Topo{1,3}.sy = sOC.Sigma_Y3(i);
    statOC.Road.Topo{1,3}.L = Ls;

    % Curviness
    % Urban
    statOC.Road.Curv{1,1}.int = sOC.Lambda_C1(i)/1000;
    statOC.Road.Curv{1,1}.mC = sOC.M_C1(i);
    statOC.Road.Curv{1,1}.sC = sOC.Sigma_C1(i);
    statOC.Road.Curv{1,1}.r_turn = 12.5;
    statOC.Road.Curv{1,1}.mL = sOC.M_L1(i);
    statOC.Road.Curv{1,1}.sL = sOC.Sigma_L1(i);
    statOC.Road.Curv{1,1}.l_turn = 2;
    % Rural
    statOC.Road.Curv{1,2}.int = sOC.Lambda_C2(i)/1000;
    statOC.Road.Curv{1,2}.mC = sOC.M_C2(i);
    statOC.Road.Curv{1,2}.sC = sOC.Sigma_C2(i);
    statOC.Road.Curv{1,2}.r_turn = 12.5;
    statOC.Road.Curv{1,2}.mL = sOC.M_L2(i);
    statOC.Road.Curv{1,2}.sL = sOC.Sigma_L2(i);
    statOC.Road.Curv{1,2}.l_turn = 2;
    % Rural
    statOC.Road.Curv{1,3}.int = sOC.Lambda_C3(i)/1000;
    statOC.Road.Curv{1,3}.mC = sOC.M_C3(i);
    statOC.Road.Curv{1,3}.sC = sOC.Sigma_C3(i);
    statOC.Road.Curv{1,3}.r_turn = 12.5;
    statOC.Road.Curv{1,3}.mL = sOC.M_L3(i);
    statOC.Road.Curv{1,3}.sL = sOC.Sigma_L3(i);
    statOC.Road.Curv{1,3}.l_turn = 2;

    % Stop signs
    % Urban
    statOC.Road.Stop{1,1}.int = sOC.Lambda_s1(i)/1000;
    statOC.Road.Stop{1,1}.t = [sOC.T_min1(i), sOC.T_min1(i)+sOC.Delta_T1(i)];
    % Rural
    statOC.Road.Stop{1,2}.int = sOC.Lambda_s2(i)/1000;
    statOC.Road.Stop{1,2}.t = [sOC.T_min2(i), sOC.T_min2(i)+sOC.Delta_T2(i)];
    % Highway
    statOC.Road.Stop{1,3}.int = sOC.Lambda_s3(i)/1000;
    statOC.Road.Stop{1,3}.t = [sOC.T_min3(i), sOC.T_min3(i)+sOC.Delta_T3(i)];

    % Speed bumps
    % Urban
    statOC.Road.Bump{1,1}.int = 0;
    statOC.Road.Bump{1,1}.v = [0, 0];
    % Rural
    statOC.Road.Bump{1,2}.int = 0;
    statOC.Road.Bump{1,2}.v = [0, 0];
    % Highway
    statOC.Road.Bump{1,3}.int = 0;
    statOC.Road.Bump{1,3}.v = [0, 0];

    % Speed signs
    % Urban
    statOC.Road.Speed{1,1}.Speed = [30, 40, 50]/3.6;
    statOC.Road.Speed{1,1}.PP = [0 sOC.p_V12_R1(i) sOC.p_V13_R1(i);
                                 sOC.p_V21_R1(i) 0 sOC.p_V23_R1(i);
                                 sOC.p_V31_R1(i) sOC.p_V32_R1(i) 0];
    statOC.Road.Speed{1,1}.LL = [sOC.L_V1_R1(i), sOC.L_V2_R1(i), sOC.L_V3_R1(i)];
    % Rural
    statOC.Road.Speed{1,2}.Speed = [60, 70, 80]/3.6;
    statOC.Road.Speed{1,2}.PP = [0 sOC.p_V12_R2(i) sOC.p_V13_R2(i);
                                 sOC.p_V21_R2(i) 0 sOC.p_V23_R2(i);
                                 sOC.p_V31_R2(i) sOC.p_V32_R2(i) 0];
    statOC.Road.Speed{1,2}.LL = [sOC.L_V1_R2(i), sOC.L_V2_R2(i), sOC.L_V3_R2(i)];
    % Highway
    statOC.Road.Speed{1,3}.Speed = [90, 100]/3.6;
    statOC.Road.Speed{1,3}.PP = [0 1;
                                 1 0];
    statOC.Road.Speed{1,3}.LL = [sOC.L_V1_R3(i), sOC.L_V2_R3(i)];

    numberToAppend = num2str(i);
    save([designationDir '\sOC_' numberToAppend], 'statOC')
end


%% ASSEMBLE MEAN sOC

designationDir = 'C:\Users\luigir\Desktop\sOC for TRD\mean_sOC';
statOC = [];

%  Road types
    statOC.Road.TC.RoadTypes = {'City','Rural','Highway'};
    statOC.Road.TC.L = Ls;
    statOC.Road.TC.LL = [1.38 6.77 7.57]*1000;
    statOC.Road.TC.PP = [0 0.89 0.11;
                         0.55 0 0.45;
                         0.11 0.89 0];
    statOC.Road.TC.Ltot = 81.72*1000;
    statOC.Road.TC.Speed = [50, 80, 100]/3.6;


    % Topography
    % Urban
    statOC.Road.Topo{1,1}.type = 'G1';
    statOC.Road.Topo{1,1}.a = sin(pi/2-2*Ls./291);
    statOC.Road.Topo{1,1}.sy = 2.30;
    statOC.Road.Topo{1,1}.L = Ls;
    % Rural
    statOC.Road.Topo{1,2}.type = 'G1';
    statOC.Road.Topo{1,2}.a = sin(pi/2-2*Ls./811);
    statOC.Road.Topo{1,2}.sy = 1.76;
    statOC.Road.Topo{1,2}.L = Ls;
    % Highway
    statOC.Road.Topo{1,3}.type = 'G1';
    statOC.Road.Topo{1,3}.a = sin(pi/2-2*Ls./739);
    statOC.Road.Topo{1,3}.sy = 1.58;
    statOC.Road.Topo{1,3}.L = Ls;

    % Curviness
    % Urban
    statOC.Road.Curv{1,1}.int = 6.59/1000;
    statOC.Road.Curv{1,1}.mC = 4.27;
    statOC.Road.Curv{1,1}.sC = 1.03;
    statOC.Road.Curv{1,1}.r_turn = 12.5;
    statOC.Road.Curv{1,1}.mL = 3.92;
    statOC.Road.Curv{1,1}.sL = 0.65;
    statOC.Road.Curv{1,1}.l_turn = 2;
    % Rural
    statOC.Road.Curv{1,2}.int = 1.92/1000;
    statOC.Road.Curv{1,2}.mC = 4.86;
    statOC.Road.Curv{1,2}.sC = 0.92;
    statOC.Road.Curv{1,2}.r_turn = 12.5;
    statOC.Road.Curv{1,2}.mL = 3.91;
    statOC.Road.Curv{1,2}.sL = 0.59;
    statOC.Road.Curv{1,2}.l_turn = 2;
    % Rural
    statOC.Road.Curv{1,3}.int = 0.89/1000;
    statOC.Road.Curv{1,3}.mC = 5.07;
    statOC.Road.Curv{1,3}.sC = 0.82;
    statOC.Road.Curv{1,3}.r_turn = 12.5;
    statOC.Road.Curv{1,3}.mL = 3.79;
    statOC.Road.Curv{1,3}.sL = 0.54;
    statOC.Road.Curv{1,3}.l_turn = 2;

    % Stop signs
    % Urban
    statOC.Road.Stop{1,1}.int = 1.18/1000;
    statOC.Road.Stop{1,1}.t = [21, 277];
    % Rural
    statOC.Road.Stop{1,2}.int = 0.1/1000;
    statOC.Road.Stop{1,2}.t = [19.75, 95.19];
    % Highway
    statOC.Road.Stop{1,3}.int = 0.06/1000;
    statOC.Road.Stop{1,3}.t = [34.50, 64.79];

    % Speed bumps
    % Urban
    statOC.Road.Bump{1,1}.int = 0;
    statOC.Road.Bump{1,1}.v = [0, 0];
    % Rural
    statOC.Road.Bump{1,2}.int = 0;
    statOC.Road.Bump{1,2}.v = [0, 0];
    % Highway
    statOC.Road.Bump{1,3}.int = 0;
    statOC.Road.Bump{1,3}.v = [0, 0];

    % Speed signs
    % Urban
    statOC.Road.Speed{1,1}.Speed = [30, 40, 50]/3.6;
    statOC.Road.Speed{1,1}.PP = [0 0.48 0.52;
                                 0.49 0 0.51;
                                 0.45 0.55 0];
    statOC.Road.Speed{1,1}.LL = [0.54 0.49 1.44]*1000;
    % Rural
    statOC.Road.Speed{1,2}.Speed = [60, 70, 80]/3.6;
    statOC.Road.Speed{1,2}.PP = [0 0.54 0.46;
                                 0.32 0 0.68;
                                 0.22 0.78 0];
    statOC.Road.Speed{1,2}.LL = [0.97 1.36 6.47]*1000;
    % Highway
    statOC.Road.Speed{1,3}.Speed = [90, 100]/3.6;
    statOC.Road.Speed{1,3}.PP = [0 1;
                                 1 0];
    statOC.Road.Speed{1,3}.LL = [7.0 1.83]*1000;

    save([designationDir '\sOC'], 'statOC')