% Program to determine at what time of day the fast charge starts at

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inputs
% Input the EV range in miles
Range = 200;

% Input the terminus park duration in format: (hour,minutes,seconds)
Terminus_Park = duration(8,0,0);

% Two dates to analyse between, if the entire survey is to be used start at
% 2000 and end date 2010 will work as inputs
% First date
Date1 = '23/11/2000';
% Second date
Date2 = '29/11/2020';

% Input 0 for fast charges by time of day and 1 for day of week
InputTD = 0;

% If the input above is 1 for day of the week the hour distribution can be
% found for the day of the week specified here, Sunday is 1, Saturday is 7,
% if 0 is entered a plot of day of week versus distribution will be
% displayed
InputD = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create vector for the journey duration
JourneyDuration = (StartEnd(:,2) - StartEnd(:,1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create a vector for the parking duration

% Create a temporary vector for the park duration, equal to the start time
% of the journey minus the finish time of the previous journey, to subtract
% the previous journey from the current journey the start and end vectors
% need to be shifted, one data point is therefore added to the end of the
% start time and one data point added to the beginning of the end time
% vector
ParkDuration = [StartEnd(:,1);datetime(datevec(1))] - [datetime(datevec(1));StartEnd(:,2)];

% The first and last points of the vector are therefore not relevant,
% remove the last data point and replace the first datapoint with NaN
ParkDuration = ParkDuration(1:length(IdentifierDistanceEndloc));
ParkDuration(1) = NaN;

% A vector is created of zeros and ones, a one is placed at the first data
% point of a new car
FirstPoint = any(circshift(diff(IdentifierDistanceEndloc(:,1)),1),length(IdentifierDistanceEndloc));
FirstPoint(1) = 1;
FirstPoint(length(IdentifierDistanceEndloc))=0;

% At all the first points put a NaN instead of the park duration
ParkDuration(FirstPoint==1)=NaN;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine time that fast charge starts

% Create a vector for the fast charge times
FC = NaT(length(IdentifierDistanceEndloc),1);

% Put the distance in the journey as the tempsum value
tempsum = IdentifierDistanceEndloc(1,2);

% Set a value for x, which increases by 1 everytime a fast charge occurs
NumFastCharges = 0;

% Go through all values in the matrix
for i = 2:length(IdentifierDistanceEndloc)
    
    % Set a count value, this value is 0 unless a fast charge is required
    % twice for one journey
    count = 0;
    
    % If the park duration is less than the terminus park duration then add
    % the distance to tempsum value
    if IdentifierDistanceEndloc(i-1,3) < 0.5 && ParkDuration(i) < Terminus_Park
        
        tempsum = tempsum + IdentifierDistanceEndloc(i,2);
    
    % Else tempsum is equal to the distance of the journey
    else
        tempsum = IdentifierDistanceEndloc(i,2);
    
    end
    
    % If the tempsum is then greater than the specified range a fast
    % charge is required
    while tempsum > Range
                
        % Calculate the average speed in miles per hour
        tempSpeed = IdentifierDistanceEndloc(i,2) / hours(JourneyDuration(i));
        
        % If count is 0 it is the first time a fast charge is conducted on
        % the journey, if count is greater than 0 it means that a fast
        % charge has already been conducted on the journey
        if count == 0
            
            % Calculate after how many miles the fast charge occurs
            tempDistance1 = Range - (tempsum - IdentifierDistanceEndloc(i,2));
            
            % Calculate the time to the fast charge in hours
            tempTimeToFC = tempDistance1 / tempSpeed;

            % Add the time to the start of the journey
            FC(NumFastCharges+1) = StartEnd(i,1) + hours(tempTimeToFC);
            
            % Set the tempsum value to the remaining distance to travel
            tempsum = IdentifierDistanceEndloc(i,2) - tempDistance1;
            
        % Else for the case if a fast charge has already occured    
        else
            
            % As this is the second fast charge of the journey it occurs
            % after the range distance, calculate the time to the fast
            % charge in hours
            tempTimeToFC = Range / tempSpeed;
            
            % Add the time to the previous fast charge time
            FC(NumFastCharges+1) = FC(NumFastCharges) + hours(tempTimeToFC);
            
            % Calculate the remaining distance to travel after the fast
            % charge
            tempsum = IdentifierDistanceEndloc(i,2) - tempDistance1 - count * Range;
            
        end
        
        % Add 1 to the number of fast charges that have occured
        NumFastCharges = NumFastCharges + 1;        
               
        % Add 1 to the count value
        count = count + 1;

    end
    
end

% Set the display format for the FC vector, if this is not done time will
% not be displayed when viewing the vector values
FC.Format = 'default';

% Reduce the size of the fast charge vector to only include elements which
% have a value
FC = FC(1:NumFastCharges);

clear count i tempDistance1 tempSpeed tempsum tempTimeToFC

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine the hour of the day that fast charges occur

% Put the dates that were read in in the correct format
Date1 = strcat(Date1,' 00:00:01');
Date2 = strcat(Date2,' 23:59:59');
Date1 = datetime(Date1,'InputFormat','dd/MM/yyyy HH:mm:ss');
Date2 = datetime(Date2,'InputFormat','dd/MM/yyyy HH:mm:ss');

% Change the Fast charging vector to only include fast charges that occur
% between the dates specified
FC = FC(FC > Date1 & FC < Date2);

% Create a vector of the hour or day or the day of the week that fast
% charge starts at depending on input preferences
if InputTD < 0.5
    FCStarts = hour(FC);
elseif InputD < 0.5
    FCStarts = weekday(FC);
else
    hourFCStarts = hour(FC);
    dayFCStarts = weekday(FC);
    FCStarts = hourFCStarts(dayFCStarts == InputD);
end
    
% Create a vector which counts the frequency with which a fast charge
% occurs in each hour, the vector will have 24 rows if finding by time of
% day or 7 rows if finding by day of week
freq = histcounts(FCStarts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine the number of days in the study

% Create 2 vectors for the start and end days of when each car joined the
% study
StartDays = StartEnd(FirstPoint,1);
EndDays = StartEnd(circshift(FirstPoint,-1),2);

% Number of cars in the study, NumCars is the entire number of cars in the
% study, while NumCarsBet is the number of cars active in the study between
% the two dates specified
NumCars = sum(FirstPoint);
NumCarsBet = sum(StartDays < Date1 & EndDays > Date2);

% Calculate the total days in the  study

% If looking between two dates this is equal to the number of cars active
% between the two dates multiplied by the number of days between the two
% dates
if Date1 > min(StartDays) || Date2 < max(EndDays)
    DaysDuration = NumCarsBet * daysact(Date1,Date2);
% Else for the case where all days in the survey are used
else
    DaysDuration = sum(days(EndDays - StartDays));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the frequency for 1 million cars
Freq1Million = ((freq / sum(DaysDuration)) * 1000000 * (InputTD * 6 + 1))';