%% start fresh
clc , clear, close all

%% Read in data
% %comment out read
Time = readmatrix('WLTP-for assignment.xlsx',... % the file
    'sheet','WLTC_class_3'... % the sheet
    ,'range','C8:C1808'); % the range

Velocity = readmatrix('WLTP-for assignment.xlsx',... % the file
    'sheet','WLTC_class_3'... % the sheet
    ,'range','E8:E1808'); % the range

Velocity = Velocity/3.6; % turns it in to m/s
%% save as mat file

%% create load
load('WLPTVelcoityTime.mat')
%% calculate acceleration
Accel = nan(size(Time)); % creates variable of NaNs
Accel(2:end,:) = diff(Velocity); % m/s as time step is one s

%% Aero force

rho = 1.202; % kg/m^3
A = 2.22; % frontal Area m^2
Cd = 0.274; % coefficent of drag per m^2

ForceAreo = 0.5.*rho.*A.*Cd.*(Velocity.^2); % N force

Fr = 0.015; % coefficient of friction
mass = 1776; % vehicle mass
ForceRR = Fr * ( mass*9.81); % 

%% Inertia loss

ForceAccel = mass*Accel;

%% Total Resistance

ForceTotal = ForceAreo+ForceRR+ForceAccel;

figure;
plot(Time,ForceAreo,...
    Time,repmat(ForceRR,size(ForceAreo)),...
    Time,ForceAccel,Time,ForceTotal)
ylabel('force N')
xlabel('Time(s)')
legend('Aero','Rolling','Accel','Total')
grid on

%% Energy per second
WheelTorque = ForceTotal.*Velocity; % joules
%% used energy
index = WheelTorque>0;
EnergyTractive = zeros(size(Time));
EnergyTractive(index) = WheelTorque(index);
%% regen energy
index = WheelTorque<0;
EnergyRegen = zeros(size(Time));
EnergyRegen(index) = WheelTorque(index);
%% total energy needed for 1 cycle excluding regen
energy1cycle = cumsum(EnergyTractive);

%% number of cycles
batteryTot = 75.*3.6e6;
NoCycles = batteryTot/energy1cycle(end,end);

%% range WLTP km
distance = cumsum(Velocity);
WLTPrangekm = NoCycles.*distance(end,end)/1000;

%% Wh/mi
WhPmi = ((energy1cycle(end,end)/3600)./(distance(end,end)/1609))/0.6;

%% energy store level 
Battery = ones(size(Accel)).*batteryTot;
for ii = 1:length(Accel)
    Battery(ii) = Battery(ii)-energy1cycle(ii);
    
    if Battery(ii)<0
        
        Battery(ii) = 0;
          
    end
end