%% Import Results GT SUITE

% The following code aims to import the results generated in GT-SUITE and 
% Open Modelica generated during a Masters Thesis entitled "Vehicle 
% Performance Simulations with OpenModelica"

%% Preliminaries

clear
close all
clc

%% Select File Title

disp("Select A Software from the Menu...")
choice1 = menu("What Software is the Data From?", ["GT SUITE", "OMedit", "Both"]);
if choice1 == 1
    s1 = "GT";
    s0 = ".txt";
elseif choice1 == 2
    s1 = "OM";
    s0 = ".csv";
else
    s1 = "Both";
    s0 = ".txt"; 
end
disp(join(["Selection: " s1]))

disp("Select the Powertrain Architecture...")
choice2 = menu("What Transmission are we Referring to?", ["BEV", "FCH"]);
if choice2 == 1
    s2 = "BEV";
    disp(join(["Selection: " s2]))
    disp(" ")
    disp("Conditions: ")
    
    choice3 = menu("What is the Simulation Goal?", ["Fuel Economy", "Acceleration/Slopes"]);
    if choice3 == 1
        s3 = "";
        choice4 = menu("What is the Simulated Driving Cycle?", ["SORT1", "SORT2", "SORT3", "Braunschweig", "MLTB", "Gillingham", "Gillingham x3"]);
        if choice4 == 1
            s4 = "SORT1";
        elseif choice4 == 2
            s4 = "SORT2";
        elseif choice4 == 3
            s4 = "SORT3";
        elseif choice4 == 4
            s4 = "Braunschweig";
        elseif choice4 == 5
            s4 = "MLTB";
        elseif choice4 == 6
            s4 = "Gillingham";
        elseif choice4 == 7
            s4 = "Gillinghamx3";
        end
        disp(join(["Selection: " s4 "(Speed Cycle)"]))
        s5 = "80";
        disp(join(["Selection: (def) " s5 "% (SOC0)"]))
        s6 = "";
        disp("Selection: (def) 7000 W (Auxiliaries)")
    else
        s3 = "Accel_";
        choice4 = menu("What is the Slope?", ["0 deg", "4 deg", "8 deg"]);
        if choice4 == 1
            s4 = "0deg";
        elseif choice4 == 2
            s4 = "4deg";
        elseif choice4 == 3
            s4 = "8deg";
        end
        disp(join(["Selection: " s4 "(Slope)"]))
        
        choice5 = menu("Inital SOC?", ["85 %", "65 %", "45 %", "35 %", "25 %"]);
        if choice5 == 1
            s5 = "85";
        elseif choice5 == 2
            s5 = "65";
        elseif choice5 == 3
            s5 = "45";
        elseif choice5 == 4
            s5 = "35";
        elseif choice5 == 5
            s5 = "25";
        end
        disp(join(["Selection: " s5 "% (SOC0)"]))
        s6 = "_Aux75";
        disp("Selection: 7500 W (Auxiliaries)")
    end
else
    s2 = "FCH";
    disp(join(["Selection: " s2]))
    disp(" ")
    disp("Conditions: ")

    choice3 = menu("What is the Simulation Goal?", ["Fuel Economy", "Acceleration/Slopes"]);
    if choice3 == 1
        s3 = "";
        choice4 = menu("What is the Simulated Driving Cycle?", ["SORT1", "SORT2", "SORT3", "Braunschweig", "MLTB", "Gillingham", "Gillingham x3"]);
        if choice4 == 1
            s4 = "SORT1";
        elseif choice4 == 2
            s4 = "SORT2";
        elseif choice4 == 3
            s4 = "SORT3";
        elseif choice4 == 4
            s4 = "Braunschweig";
        elseif choice4 == 5
            s4 = "MLTB";
        elseif choice4 == 6
            s4 = "Gillingham";
        elseif choice4 == 7
            s4 = "Gillinghamx3";
        end
        disp(join(["Selection: " s4 "(Speed Cycle)"]))

        choice5 = menu("Initial SOC?", ["62 %", "30 %"]);
        if choice5 == 1
            s5 = "62";
        else
            s5 = "30";
        end
        disp(join(["Selection: " s5 "% (SOC0)"]))
        
        choice6 = menu("Auxiliaries?", ["7k W", "0 W"]);
        if choice6 == 1
            s6 = "";
            disp("Selection: 7000 W (Auxiliaries)")
        else
            s6 = "_Aux0";
            disp("Selection: 0 W (Auxiliaries)")
        end
    else
        s3 = "Accel_";
        choice4 = menu("What is the Slope?", ["0 deg", "2 deg","4 deg", "6 deg", "8 deg", "10 deg", "12 deg"]);
        if choice4 == 1
            s4 = "0deg";
        elseif choice4 == 2
            s4 = "2deg";
        elseif choice4 == 3
            s4 = "4deg";
        elseif choice4 == 4
            s4 = "6deg";
        elseif choice4 == 5
            s4 = "8deg";
        elseif choice4 == 6
            s4 = "10deg";
        elseif choice4 == 7
            s4 = "12deg";
        end
        disp(join(["Selection: " s4 "(Slope)"]))

        choice5 = menu("Initial SOC?", ["62 %", "30 %"]);
        if choice5 == 1
            s5 = "62";
        else
            s5 = "30";
        end
        disp(strjoin(["Selection: " s5 " % (SOC0)"]))
        s6 = "_Aux75";
        disp("Selection: 7500 W (Auxiliaries)")
    end
end

choice7 = menu("Vehicle Mass?", ["15 570 kg", "19 000 kg"]);
if choice7 == 1
    s7 = "15t";
    disp("Selection: 15 570 kg")
else
    s7 = "19t";
    disp("Selection: 19 000 kg")
end

temp = strjoin([s3 s2 "_" s1 "_" s4 "_" s7 "_" s5 s6 s0]);
fileName = strrep(temp, " ", "");


%% Define Variables

if choice1 == 1
    Data1 = importdata(fileName,',', 4);
elseif choice1 == 2
    Data1 = importdata(fileName,',', 2);
else
    Data1 = importdata(strrep(strjoin([s3 s2 "_GT_" s4 "_" s7 "_" s5 s6 ".txt"]), " ", ""));
    Data2 = importdata(strrep(strjoin([s3 s2 "_OM_" s4 "_" s7 "_" s5 s6 ".csv"]), " ", ""));
end
clear temp s0 s2 s3 s4 s6 s7

time1 = Data1.data(:, 1);               % [s] Time Elapsed
RefSpeed1 = Data1.data(:, 2);           % [m/s] Reference Velocity
Inc1 = Data1.data(:, 3);                % [deg] Inclination
RealSpeed1 = Data1.data(:, 4);          % [m/s] Real Velocity
Throttle1 = Data1.data(:, 5);           % [-] Throttle Pedal
Brake1 = Data1.data(:, 6);              % [-] Brake Pedal
SOC1 = Data1.data(:, 7);                % [-] State of Charge
OCV1 = Data1.data(:, 8);                % [V] Open Circuit Voltage
Req1 = Data1.data(:, 9);                % [Ω] Equivalent Internal Resistance
Vbatt1 = Data1.data(:, 10);             % [V] Battery Terminal Voltage
Ibatt1 = Data1.data(:, 11);             % [A] Battery Current
if size(Data1.data, 2) > 11
    Iaux1 = Data1.data(:, 12);          % [A] Auxiliary Current
    Torque1 = Data1.data(:, 13);        % [Nm] Motor Torque
    RPM1 = Data1.data(:, 14);           % [RPM] Motor Rotational Velocity
    torque_flag = 0;
else
    torque_flag = 1;
    disp(' ')
    disp('### MISSING TORQUE DATA ###')
end
if choice2 == 2
    State1 = Data1.data(:,15);          % [-] Fuel Cell Control State
    H2_dot1 = Data1.data(:,16);         % [g/s] Hydrogen Consumption Rate
    H21 = Data1.data(:, 17);            % [g] Hydrogen Consumption
    Vfc1 = Data1.data(:,18);            % [V] Fuel Cell Voltage
    Ifc1 = Data1.data(:, 19);           % [A] Fuel Cell Current
end

if choice1 == 3
    time2 = Data2.data(:, 1);           % [s] Time Elapsed
    RefSpeed2 = Data2.data(:, 2);       % [m/s] Reference Velocity
    Inc2 = Data2.data(:, 3);            % [deg] Inclination
    RealSpeed2 = Data2.data(:, 4);      % [m/s] Real Velocity
    Throttle2 = Data2.data(:, 5);       % [-] Throttle Pedal
    Brake2 = Data2.data(:, 6);          % [-] Brake Pedal
    SOC2 = Data2.data(:, 7);            % [-] State of Charge
    OCV2 = Data2.data(:, 8);            % [V] Open Circuit Voltage
    Req2 = Data2.data(:, 9);            % [Ω] Equivalent Internal Resistance
    Vbatt2 = Data2.data(:, 10);         % [V] Battery Terminal Voltage
    Ibatt2 = Data2.data(:, 11);         % [A] Battery Current
    Iaux2 = Data2.data(:, 12);          % [A] Auxiliary Current
    Torque2 = Data2.data(:, 13);        % [Nm] Motor Torque
    RPM2 = Data2.data(:, 14);           % [RPM] Motor Rotational Velocity
    if choice2 == 2
        State2 = Data2.data(:,15);      % [-] Fuel Cell Control State
        H2_dot2 = Data2.data(:,16);     % [kg/s] Hydrogen Consumption Rate
        H22 = Data2.data(:, 17);        % [kg] Hydrogen Consumption
        Vfc2 = Data2.data(:,18);        % [V] Fuel Cell Voltage
        Ifc2 = Data2.data(:, 19);       % [A] Fuel Cell Current
    end
end

if choice1 == 1 || choice1 == 3
    RealSpeed1 = RealSpeed1/3.6;
    RefSpeed1 = RefSpeed1/3.6;
end

flag = 0;
% The following if condition fixes some files who had the acceleration exported instead of the velocity.
dt = time1(2, 1) - time1(1, 1);
if (choice1 == 1 || choice1 == 3) && abs(sum(RealSpeed1)/size(RealSpeed1,1)) < 1 && choice3 == 1
    flag = 1;
    disp(" ")
    disp("DATA CONTAINED ERRORS!")
    disp("⠀⠀⠀⠀⣀⠤⠔⠒⠒⠒⠒⠒⠒⠒⠦⢄⣀⠀⠀⠀⠀")
    disp("⠀⢀⡴⠋⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⢄⠀⠀")
    disp("⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢣⠀")
    disp("⢸⠀⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢢⠈⡇")
    disp("⢸⠀⣇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡼⠀⡇")
    disp("⠘⡆⢸⠀⢀⣀⣤⣄⡀⠀⠀⠀⢀⣤⣤⣄⡀⠀⡇⡸⠀")
    disp("⠀⠘⣾⠀⣿⣿⣿⣿⣿⠀⠀⠀⣿⣿⣿⣿⣿⠀⡗⠁⠀")
    disp("⠀⠀⣿⠀⠙⢿⣿⠿⠃⢠⢠⡀⠙⠿⣿⠿⠃⠀⡇⠀⠀")
    disp("⠀⠀⠘⣄⡀⠀⠀⠀⢠⣿⢸⣿⠀⠀⠀⠀⠀⣠⠇⠀⠀")
    disp("⠀⠀⠀⠀⡏⢷⡄⠀⠘⠟⠈⠿⠁⠀⢠⡞⡹⠁⠀⠀⠀")
    disp("⠀⠀⠀⠀⢹⠸⠘⢢⢠⠤⠤⡤⡄⢰⢡⠁⡇⠀⠀⠀⠀")
    disp("⠀⠀⠀⠀⢸⠀⠣⣹⢸⠒⠒⡗⡇⣩⠌⢀⡇⠀⠀⠀⠀")
    disp("⠀⠀⠀⠀⠈⢧⡀⠀⠉⠉⠉⠉⠁⠀⣀⠜⠀⠀⠀⠀⠀")
    disp("⠀⠀⠀⠀⠀⠀⠉⠓⠢⠤⠤⠤⠔⠊⠁⠀⠀⠀⠀⠀⠀")
    disp("DATA CONTAINED ERRORS!")
end
clear count1

%% Error Analysis

if choice3 == 1
    disp(" ")
    cycledist = sum(RefSpeed1.*time1)/1000;
    disp(strjoin(["Cycle Reference Distance" num2str(cycledist) " km"]))

    sz1 = size(time1,1);
    disp(" ")
    if choice1 ~= 3
        disp(strjoin(["Driver Errors in" s1]))
    else
        disp("Driver Errors in GT")
    end
    Err1a = sum(abs((time1(2:sz1) - time1(1:(sz1 - 1))).*(RefSpeed1(2:sz1)-RealSpeed1(2:sz1))));
    disp(strjoin(["Err1: " num2str(Err1a)]))
    Err2a = sum(abs((time1(2:sz1) - time1(1:(sz1 - 1))).*(RefSpeed1(2:sz1)-RealSpeed1(2:sz1))))/(time1(sz1));
    disp(strjoin(["Err2: " num2str(Err2a)]))
    Err3a = sum((time1(2:sz1) - time1(1:(sz1 - 1))).*(RefSpeed1(2:sz1)-RealSpeed1(2:sz1)));
    disp(strjoin(["Err3: " num2str(Err3a)]))

    if choice1 == 3
        sz2 = size(time2,1);
        disp(" ")
        disp("Driver Errors in OM")
        Err1b = sum(abs((time2(2:sz2) - time2(1:(sz2 - 1))).*(RefSpeed2(2:sz2)-RealSpeed2(2:sz2))));
        disp(strjoin(["Err1: " num2str(Err1b)]))
        Err2b = sum(abs((time2(2:sz2) - time2(1:(sz2 - 1))).*(RefSpeed2(2:sz2)-RealSpeed2(2:sz2))))/(time2(sz2));
        disp(strjoin(["Err2: " num2str(Err2b)]))
        Err3b = sum((time2(2:sz2) - time2(1:(sz2 - 1))).*(RefSpeed2(2:sz2)-RealSpeed2(2:sz2)));
        disp(strjoin(["Err3: " num2str(Err3b)]))
    end
else
    if choice1 == 1 %GT
        sz1 = size(time1, 1);
        tgt_vect = 50*ones(sz1, 1);
        [~, idx] = min(abs(3.6*RealSpeed1 - tgt_vect));
        ZtoFty_1 = time1(idx, 1);
        disp(" ")
        disp(['Acceleration Time (0-50 kph): ' num2str(ZtoFty_1) ' s (GT-SUITE)'])
    elseif choice1 == 2 %OM
        sz1 = size(time1, 1);
        tgt_vect = 50*ones(sz1, 1);
        [~, idx] = min(abs(3.6*RealSpeed1 - tgt_vect));
        ZtoFty_1 = time1(idx, 1);
        disp(" ")
        disp(['Acceleration Time (0-50 kph): ' num2str(ZtoFty_1) ' s (OMedit)'])
    else %Both
        sz1 = size(time1, 1);
        tgt_vect = 50*ones(sz1, 1);
        [~, idx] = min(abs(3.6*RealSpeed1 - tgt_vect));
        ZtoFty_1 = time1(idx, 1);
        disp(" ")
        disp(['Acceleration Time (0-50 kph): ' num2str(ZtoFty_1) ' s (GT-SUITE)'])
        clear tgt_vect
        sz2 = size(time2, 1);
        tgt_vect = 50*ones(sz2, 1);
        [~, idx] = min(abs(3.6*RealSpeed2 - tgt_vect));
        ZtoFty_1 = time2(idx, 1);
        disp(['Acceleration Time (0-50 kph): ' num2str(ZtoFty_1) ' s (OMedit)'])
    end
end

if flag == 1
    disp(" ")
    disp("DATA CONTAINED ERRORS & HAS BEEN ALTERED!")
    disp("⠀⠀⠀⠀⣀⠤⠔⠒⠒⠒⠒⠒⠒⠒⠦⢄⣀⠀⠀⠀⠀")
    disp("⠀⢀⡴⠋⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⢄⠀⠀")
    disp("⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢣⠀")
    disp("⢸⠀⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢢⠈⡇")
    disp("⢸⠀⣇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡼⠀⡇")
    disp("⠘⡆⢸⠀⢀⣀⣤⣄⡀⠀⠀⠀⢀⣤⣤⣄⡀⠀⡇⡸⠀")
    disp("⠀⠘⣾⠀⣿⣿⣿⣿⣿⠀⠀⠀⣿⣿⣿⣿⣿⠀⡗⠁⠀")
    disp("⠀⠀⣿⠀⠙⢿⣿⠿⠃⢠⢠⡀⠙⠿⣿⠿⠃⠀⡇⠀⠀")
    disp("⠀⠀⠘⣄⡀⠀⠀⠀⢠⣿⢸⣿⠀⠀⠀⠀⠀⣠⠇⠀⠀")
    disp("⠀⠀⠀⠀⡏⢷⡄⠀⠘⠟⠈⠿⠁⠀⢠⡞⡹⠁⠀⠀⠀")
    disp("⠀⠀⠀⠀⢹⠸⠘⢢⢠⠤⠤⡤⡄⢰⢡⠁⡇⠀⠀⠀⠀")
    disp("⠀⠀⠀⠀⢸⠀⠣⣹⢸⠒⠒⡗⡇⣩⠌⢀⡇⠀⠀⠀⠀")
    disp("⠀⠀⠀⠀⠈⢧⡀⠀⠉⠉⠉⠉⠁⠀⣀⠜⠀⠀⠀⠀⠀")
    disp("⠀⠀⠀⠀⠀⠀⠉⠓⠢⠤⠤⠤⠔⠊⠁⠀⠀⠀⠀⠀⠀")
end

disp(' ')
SOC0 = str2num(s5)/100;
if choice2 == 1
    if choice1 == 1
        BattEnergy1 = (SOC0 - SOC1(sz1))*475.4;
        disp(['GT Battery Energy: ' num2str(BattEnergy1) ' kWh'])
    elseif choice1 == 2
        BattEnergy1 = (SOC0 - SOC1(sz1))*475.4;
        disp(['OM Battery Energy: ' num2str(BattEnergy1) ' kWh'])
    else
        BattEnergy1 = (SOC0 - SOC1(sz1))*475.4;
        disp(['GT Battery Energy: ' num2str(BattEnergy1) ' kWh'])
        BattEnergy2 = (SOC0 - SOC2(sz2))*475.4;
        disp(['OM Battery Energy: ' num2str(BattEnergy2) ' kWh'])
        pervar1 = ((BattEnergy2 - BattEnergy1)/BattEnergy1)*100;
        disp(' ')
        disp(['Battery Energy Variation: ' num2str(pervar1) ' %'])
    end
else
    if choice1 == 1
        BattEnergy1 = (SOC0 - SOC1(sz1))*39.6;
        disp(['GT Battery Energy: ' num2str(BattEnergy1) ' kWh'])
        H21f = H21(sz1);
        disp(['GT Hydrogen Consumed: ' num2str(H21f) ' kg' ])
    elseif choice1 == 2
        BattEnergy1 = (SOC0 - SOC1(sz1))*39.6;
        disp(['OM Battery Energy: ' num2str(BattEnergy1) ' kWh'])
        H21f = H21(sz1)*1000;
        disp(['OM Hydrogen Consumed: ' num2str(H21f) ' kg' ])
    else
        BattEnergy1 = (SOC0 - SOC1(sz1))*39.6;
        disp(['GT Battery Energy: ' num2str(BattEnergy1) ' kWh'])
        H21f = H21(sz1);
        disp(['GT Hydrogen Consumed: ' num2str(H21f) ' g' ])
        BattEnergy2 = (SOC0 - SOC2(sz2))*39.6;
        disp(['OM Battery Energy: ' num2str(BattEnergy2) ' kWh'])
        H22f = H22(sz2)*1000;
        disp(['OM Hydrogen Consumed: ' num2str(H22f) ' g'])
        disp(' ')
        pervar1 = ((BattEnergy2 - BattEnergy1)/BattEnergy1)*100;
        disp(['Battery Energy Variation: ' num2str(pervar1) ' %'])
        pervar2 = ((H22f - H21f)/H21f)*100;
        disp(['Hydrogen Variation: ' num2str(pervar2) ' %'])
        EqEnergy1 = BattEnergy1 + (0.98*0.5*119000*H21f)/(1000*3600);
        disp(' ')
        disp(['GT Equivalent Energy Consumption: ' num2str(EqEnergy1) ' kWh'])
        EqEnergy2 = BattEnergy2 + (0.98*0.5*119000*H22f)/(1000*3600);
        disp(['OM Equivalent Energy Consumption: ' num2str(EqEnergy2) ' kWh'])
        disp(' ')
        pervar3 = ((EqEnergy2 - EqEnergy1)/EqEnergy1)*100;
        disp(['Equivalent Energy Variation: ' num2str(pervar3) ' %'])
    end
end

%% Fuel Economy Graphs

if choice3 == 1
% Speed Plot
figure
hold on;
grid on;
if choice1 ~= 3
    if choice1 == 1
        plot(time1, 3.6*RefSpeed1, 'k')
        plot(time1, 3.6*RefSpeed1 + 2, 'Color', "#CCCCCC")
        plot(time1, 3.6*RefSpeed1 - 2, 'Color', "#CCCCCC")
        plot(time1, 3.6*RealSpeed1, 'r--')
        set(gca, 'fontsize', 14)
    else
        plot(time1, 3.6*RefSpeed1, 'k')
        plot(time1, 3.6*RefSpeed1 + 2, 'Color', "#CCCCCC")
        plot(time1, 3.6*RefSpeed1 - 2, 'Color', "#CCCCCC")
        plot(time1, 3.6*RealSpeed1, 'b--')
        set(gca, 'fontsize', 14)
    end
    title(strjoin([s1 "Vehicle Speed Profile"]), 'FontSize', 22)
    legend('Reference Speed', '', '', 'Real Speed', 'FontSize', 16)
else 
    if flag == 0
        plot(time1, 3.6*RefSpeed1 + 2, 'Color', "#CCCCCC", 'LineWidth', 1.5)
        plot(time1, 3.6*RefSpeed1 - 2, 'Color', "#CCCCCC", 'LineWidth', 1.5)
        plot(time1, 3.6*RefSpeed1, 'k', 'LineWidth', 1.25)
        plot(time1, 3.6*RealSpeed1, 'r-.', 'LineWidth', 1.25)
        plot(time2, 3.6*RealSpeed2, '-.','LineWidth', 1.25, 'Color', [0 0.5 1])
        set(gca, 'fontsize', 14)
        title("Vehicle Speed Profile", 'FontSize', 22)
        legend('', '', 'Reference Speed', 'GT Real Speed', 'OM Real Speed', 'FontSize', 16)
    else
        plot(time1, 3.6*RefSpeed1 + 2, 'Color', "#CCCCCC", 'LineWidth', 1.5)
        plot(time1, 3.6*RefSpeed1 - 2, 'Color', "#CCCCCC", 'LineWidth', 1.5)
        plot(time1, 3.6*RefSpeed1, 'k', 'LineWidth', 1.25)
        plot(time2, 3.6*RealSpeed2, '-.','LineWidth', 1.25, 'Color', [0 0.5 1])
        set(gca, 'fontsize', 14)
        title("Vehicle Speed Profile", 'FontSize', 22)
        legend('', '', 'Reference Speed',  'OM Real Speed', 'FontSize', 16) 
    end
end
xlim([0, max(time1)])
ylim([0, 65]);
xlabel('Time [s]', 'FontSize', 16)
ylabel('Speed [km/h]', 'FontSize', 16)


% Battery Power
figure
hold on;
grid on;
if choice1 ~= 3
    if choice1 == 1
        plot(time1, Vbatt1.*Ibatt1/1000, 'r', 'LineWidth', 0.75)
        set(gca, 'fontsize', 14)
    else
        plot(time1, Vbatt1.*Ibatt1/1000, '-.','LineWidth', 0.75, 'Color', [0 0.5 1])
        set(gca, 'fontsize', 14)
    end
    title(strjoin([s1 "Battery State of Charge"]), 'FontSize', 22)
else 
    plot(time1, Vbatt1.*Ibatt1/1000, 'r-', 'LineWidth', 1.25)
    plot(time2, Vbatt2.*Ibatt2/1000, ':','LineWidth', 0.75, 'Color', [0 0.5 1])
    set(gca, 'fontsize', 14)
    title("Battery Power", 'FontSize', 22)
    legend('GT', 'OM', 'FontSize', 16)
end
xlim([0, max(time1)])
ylim([-250, 300]);
xlabel('Time [s]', 'FontSize', 16)
ylabel('Battery Power [kW]', 'FontSize', 16)

% SOC Plot
figure
hold on;
grid on;
if choice1 ~= 3
    if choice1 == 1
        plot(time1, SOC1*100, 'r', 'LineWidth', 0.75)
        set(gca, 'fontsize', 14)
    else
        plot(time1, SOC1*100, '-','LineWidth', 0.75, 'Color', [0 0.5 1])
        set(gca, 'fontsize', 14)
    end
    title(strjoin([s1 "Battery State of Charge"]), 'FontSize', 22)
else 
    plot(time1, SOC1*100, 'r-', 'LineWidth', 1.25)
    plot(time2, SOC2*100, '-','LineWidth', 1.25, 'Color', [0 0.5 1])
    set(gca, 'fontsize', 14)
    title("Battery State of Charge", 'FontSize', 22)
    legend('GT', 'OM', 'FontSize', 16)
end
xlim([0, max(time1)])
ylim([min(SOC1)*100 - 1, max(SOC1)*100 + 0.5]);
xlabel('Time [s]', 'FontSize', 16)
ylabel('SOC [%]', 'FontSize', 16)

% Torque Request Plot
figure
hold on;
grid on;
if choice1 ~= 3
    if choice1 == 1
        plot(time1, Torque1, 'r', 'LineWidth', 0.75)
        set(gca, 'fontsize', 14)
    else
        plot(time1, Torque1, '-','LineWidth', 0.75, 'Color', [0 0.5 1])
        set(gca, 'fontsize', 14)
    end
    title(strjoin([s1 "Motor Torque"]), 'FontSize', 22)
else 
        plot(time1, Torque1, 'r-', 'LineWidth', 1.25)
        plot(time2, Torque2, ':','LineWidth', 0.75, 'Color', [0 0.5 1])
        set(gca, 'fontsize', 14)
        title("Motor Torque Trace", 'FontSize', 22)
        legend('GT', 'OM', 'FontSize', 16)
end
xlim([0, max(time1)])
ylim([-1500, 1500]);
xlabel('Time [s]', 'FontSize', 16)
ylabel('Motor Torque [Nm]', 'FontSize', 16)

% Fuel Cell Operating Mode
if choice2 == 2
    figure
    hold on;
    grid on;
    if choice1 ~= 3
        if choice1 == 1
            plot(time1, State1, 'r', 'LineWidth', 0.75)
            set(gca, 'fontsize', 14)
        else
            plot(time1, State1, '-','LineWidth', 0.75, 'Color', [0 0.5 1])
            set(gca, 'fontsize', 14)
        end
        title(strjoin([s1 "Fuel Cell Operating Mode"]), 'FontSize', 22)
    else 
        plot(time2, State2, '-','LineWidth', 1.25, 'Color', [0 0.5 1])
        plot(time1, State1, 'r-.', 'LineWidth', 1.25)
        set(gca, 'fontsize', 14)
        title("Fuel Cell Operating Mode", 'FontSize', 22)
        legend('OM', 'GT', 'FontSize', 16)
    end
    xlim([0, max(time1)])
    ylim([0, 5])
    xlabel('Time [s]', 'FontSize', 16)
    ylabel('Fuel Cell Mode', 'FontSize', 16)
end

% Fuel Cell Hydrogen Consumed
if choice2 == 2
    figure
    hold on;
    grid on;
    if choice1 ~= 3
        if choice1 == 1
            plot(time1, H21, 'r', 'LineWidth', 0.75)
            set(gca, 'fontsize', 14)
            ylim([0, H21(sz1)+50])
        else
            plot(time1, H21*1000, '-','LineWidth', 0.75, 'Color', [0 0.5 1])
            set(gca, 'fontsize', 14)
            ylim([0, H21(sz1)*1000+50])
        end
        title(strjoin([s1 "Fuel Cell Hydrogen Consumption"]), 'FontSize', 22)
    else 
        plot(time1, H21, 'r-', 'LineWidth', 1.25)
        plot(time2, H22*1000, '-','LineWidth', 1.25, 'Color', [0 0.5 1])
        set(gca, 'fontsize', 14)
        title("Fuel Cell Hydrogen Consumption", 'FontSize', 22)
        legend('GT', 'OM', 'FontSize', 16)
        ylim([0, H22(sz2)*1000+50])
    end
    xlim([0, max(time1)])
    xlabel('Time [s]', 'FontSize', 16)
    ylabel('Hydrogen Consumed [g]', 'FontSize', 16)
end

% Fuel Cell Power
if choice2 == 2
    figure
    hold on;
    grid on;
    if choice1 ~= 3
        if choice1 == 1
            plot(time1, Vfc1.*Ifc1/1000, 'r', 'LineWidth', 0.75)
            set(gca, 'fontsize', 14)
        else
            plot(time1, Vfc1.*Ifc1/1000, 'b', 'LineWidth', 0.75)
            set(gca, 'fontsize', 14)
        end
        title(strjoin([s1 "Fuel Cell Hydrogen Consumption"]), 'FontSize', 22)
    else 
        plot(time1, (Vfc1.*Ifc1)/1000, 'r-', 'LineWidth', 1)
        plot(time2, (Vfc2.*Ifc2)/1000, '-.','LineWidth', 1, 'Color', [0 0.5 1])
        set(gca, 'fontsize', 14)
        title("Fuel Cell Gross Power", 'FontSize', 22)
        legend('GT', 'OM', 'FontSize', 16)
    end
    xlim([0, max(time1)])
    ylim([0, 100])
    xlabel('Time [s]', 'FontSize', 16)
    ylabel('Power [kW]', 'FontSize', 16)
end


%% Acceleration Graphs
else

figure
hold on;
grid on;
if choice1 ~= 3
    if choice1 == 1
        plot(time1, 3.6*RefSpeed1, 'k-.', 'LineWidth', 1.25)
        plot(time1, 3.6*RealSpeed1, 'r--', 'LineWidth', 1.25)
        set(gca, 'fontsize', 14)
    else
        plot(time1, 3.6*RefSpeed1, 'k-.', 'LineWidth', 1.25)
        plot(time1, 3.6*RealSpeed1, '--','LineWidth', 1.25, 'Color', [0 0.5 1])
        set(gca, 'fontsize', 14)
    end
    title(strjoin([s1 "Vehicle Speed Profile"]), 'FontSize', 22)
    legend('Reference Speed', 'Real Speed', 'FontSize', 16)
else 
        plot(time2, 3.6*RefSpeed2, 'k', 'LineWidth', 0.75)
        plot(time1, 3.6*RealSpeed1, 'r-.', 'LineWidth', 1.25)
        plot(time2, 3.6*RealSpeed2, '-.','LineWidth', 1.25, 'Color', [0 0.5 1])
        set(gca, 'fontsize', 14)
        title("Vehicle Speed Profile", 'FontSize', 22)
        legend('Reference Speed', 'GT Real Speed', 'OM Real Speed', 'FontSize', 16)
end
xlim([0, 35])%0.5*(max(time1)+max(time2))])
ylim([0, 55]);
xlabel('Time [s]', 'FontSize', 16)
ylabel('Speed [km/h]', 'FontSize', 16)


figure
hold on;
grid on;

if torque_flag == 1
    if choice1 == 1
        disp(' ')
        disp('### ERROR DISPLAYING GT TORQUE GRAPHS ###')
    else
        plot([0, 100], [1467.9 1467.9], 'k-.', 'LineWidth', 0.75)
        plot(time2, Torque2, '-.','LineWidth', 0.75, 'Color', [0 0.5 1])
        set(gca, 'fontsize', 14)
        title("Motor Torque Trace", 'FontSize', 22)
        legend('Max Torque', 'OM Motor Torque', 'FontSize', 16)
    end
else
    if choice1 ~= 3
        if choice1 == 1
            plot([0, 100], [1467.9 1467.9], 'k:')
            plot(time1, Torque1, 'r--', 'LineWidth', 1.25)
            set(gca, 'fontsize', 14)
            title("GT Motor Torque Trace", 'FontSize', 22)
        else
            plot([0, 100], [1467.9 1467.9], 'k-.')
            plot(time1, Torque1, '--','LineWidth', 1.25, 'Color', [0 0.5 1])
            set(gca, 'fontsize', 14)
            title("OM Motor Torque Trace", 'FontSize', 22)
        end
        legend('Max Torque', 'Motor Torque', 'FontSize', 16)
    else 
            plot([0, 100], [1467.9 1467.9], 'k-.', 'LineWidth', 1)
            plot(time1, Torque1, 'r-.', 'LineWidth', 1.25)
            plot(time2, Torque2, '-.','LineWidth', 1.25, 'Color', [0 0.5 1])
            set(gca, 'fontsize', 14)
            title("Motor Torque Trace", 'FontSize', 22)
            legend('Max Torque', 'GT Motor Torque', 'OM Motor Torque', 'FontSize', 16)
    end
end
xlim([0, 35])%0.5*(max(time1)+max(time2))])
ylim([0, 1500]);
xlabel('Time [s]', 'FontSize', 16)
ylabel('Torque [Nm]', 'FontSize', 16)

    if choice2 == 2
    figure
    hold on;
    grid on;
        if choice1 ~= 3
            if choice1 == 1
                plot(time1, Vfc1.*Ifc1/1000, 'r', 'LineWidth', 0.75)
                set(gca, 'fontsize', 14)
            else
                plot(time1, Vfc1.*Ifc1/1000, '-','LineWidth', 0.75, 'Color', [0 0.5 1])
                set(gca, 'fontsize', 14)
            end
            title(strjoin([s1 "Fuel Cell Hydrogen Consumption"]))
        else 
            plot(time1, (Vfc1.*Ifc1)/1000, 'r-', 'LineWidth', 1)
            plot(time2, (Vfc2.*Ifc2)/1000, '-.','LineWidth', 0.75, 'Color', [0 0.5 1])
            set(gca, 'fontsize', 14)
            title("Fuel Cell Gross Power", 'FontSize', 22)
            legend('GT', 'OM', 'FontSize', 16)
        end
        xlim([0, 35])
        ylim([0, 100])
        xlabel('Time [s]', 'FontSize', 16)
        ylabel('Power [kW]', 'FontSize', 16)
    end

    if choice2 == 2
    figure
    hold on;
    grid on;
        if choice1 ~= 3
            if choice1 == 1
                plot(time1, State1, 'r', 'LineWidth', 0.75)
                set(gca, 'fontsize', 14)
            else
                plot(time1, State1, '-','LineWidth', 0.75, 'Color', [0 0.5 1])
                set(gca, 'fontsize', 14)
            end
            title(strjoin([s1 "Fuel Cell Operating Mode"]))
        else 
            plot(time2, State2, '-','LineWidth', 1.25, 'Color', [0 0.5 1])
            plot(time1, State1, 'r-.', 'LineWidth', 1.25)
            set(gca, 'fontsize', 14)
            title("Fuel Cell Operating Mode", 'FontSize', 22)
            legend('OM', 'GT', 'FontSize', 16)
        end
        xlim([0, 35])
        ylim([0, 5])
        xlabel('Time [s]', 'FontSize', 16)
        ylabel('Fuel Cell Mode', 'FontSize', 16)
    end
end
