function [dayLog,thisDay,kStep] = doMLEPSim(Tset,dayLog,ep,TIMESTEP,dayFlag,kStep,prevkStep,logData,iter,coolDownCount,COOLDOWNDAYS,Tbase,numMultiDays,myDefaultScheduleTset)
%% Create an mlepProcess instance and start EP. (First run only).

global timeSim
global timeRestart
%======================================================================
% OPTIONS
%======================================================================
BUFFERLIMIT = 7; % 3 or 7              % DISABLE IF USING CHECKPOINT. How many days to buffer with default schedule? 
USECHECKPOINT = 0;             % Use checkpoint?
WEATHERNAME = 'USA_HI_Honolulu.911820_TMY2';
IDFNAME = 'SmallOffice2004TsetTest_V830';
NUMCOLUMNS = 3+7;          % Should be number returned parameters. Should be length(SP) + num returned parameters
NUMDAYS = 9999;           % must match number of times to repeat in IDF
TAKEMEANTIME = 1;

VERNUMBER = 2;  % version For:number of communication protocol (2 for E+ 7.2.0)
MODE = 1;       % 0 = COMPLETELY OFF.  1 = Normal Tset mode.  2 = Duty Cycle mode (do not use).
DOCOAST = 1;    % Do coasting behavior if Tset > Tamb

%======================================================================


NUMSTEPS = 60/TIMESTEP;  % Matches numsteps/hr as stated in IDF exactly
DELTATIME = TIMESTEP*60;   
MAXSTEPS = NUMSTEPS*24*NUMDAYS*(BUFFERLIMIT+1);  
DEADZONE = 0.1;
BUFFERDAY = 0;                 % Incrementer to track buffered days.
HVAC_ON = 1; % For first hour
% Schedules to check:
%------------------------------------
% CoolingCoilAvailSched
% Min OA Sched
% These passively contribute to HVAC power consumption. These schedules may
% need to be changed, since the water must be cooled before it is availble
% for use. For now:

% COOLING OFF: 12AM - 6AM; 8PM - 12AM
% COOLING ON: 6AM - 8PM

%% The main simulation loop
% MAXSTEPS to have +1 to make the vectors line up, but we can disregard the
% last sample. The kStep marking each day should be (MAXSTEPS-1)/NUMDAYS
    id = get(getCurrentTask,'ID');
    if isempty(id)  % If ID is still empty after the above assignment, you're running serial. 
        id = 1;     % Make sure we don't have an empty index.
    end 

    if iter(id)  == 1 && coolDownCount == 1
	timeSim = zeros(300,length(iter)); % Store time averages.
	timeRestart = timeSim;
    end
    
    % Below: the amount of buffer days is set to be the limit to emulate
    % the state of all buffer days being done. It is conditionally set by
    % two important conditions:
    % 1) The number of cooldown days is positive; these days are when the default
    %    Tset schedule is used to incooporate variable terminal costs and
    %    to achieve a steady state condition. 
    % 2) The number of iterations exceeds 1. Because cooldown days help in
    %    achieving steady state, there is no need to use redundant buffer days.
    if coolDownCount > 0 || iter(id) > 1, BUFFERDAY = BUFFERLIMIT; end 
    %fprintf('BufferDay = %d / %d\n',BUFFERDAY,BUFFERLIMIT)
    fprintf('coolDownCount = %d / %d ',coolDownCount,COOLDOWNDAYS)
    if coolDownCount < numMultiDays
        fprintf('(in planning horizon)')
    else
        fprintf('(in termination horizon)')
    end
    fprintf('\n')
    
 
TsetTest = Tset;               % Tset we aim to test with this iteration
%Tset = myDefaultScheduleTset;  % Existing schedule (previous day's) to buffer with.  

while kStep <= MAXSTEPS    

    if BUFFERDAY == BUFFERLIMIT                         % Buffer with default schedule until
        Tset = TsetTest;                                % we're done buffering. Then record
    else                                                 % with ACTUAL Tset. 
        Tset = myDefaultScheduleTset(BUFFERDAY+1,:);
    end
    % TsetTest takes on the default schedule as specified in mainEP. 

    packet = ep.read;       % Read a data packet from EP.
    if isempty(packet), error('Could not read outputs from E+.'); end
    [mlepflag, eptime, outputs] = mlepDecodePacket(packet);     % Parse it to obtain building outputs
    if mlepflag ~= 0, break; end % Error in simulation or the simulation already ended.
    if USECHECKPOINT == 1    
        if kStep == 1 && iter(id) == 1 && coolDownCount == 0
            disp(['Initial checkpoint, kStep=' num2str(kStep) ' MAXSTEPS=' num2str(MAXSTEPS), ' Worker=' num2str(id)])
            ep.checkpoint()
            disp('First successful restart.')
            tic
        elseif kStep == 1 && iter(id) > 1 && coolDownCount == 0
            disp(['Restart, kStep=' num2str(kStep+288*COOLDOWNDAYS) ' MAXSTEPS=' num2str(MAXSTEPS), ' Worker=' num2str(id)])
            timeSim(iter(id),id) = toc;
            fprintf(['Simulation Time Elapsed: ', num2str(timeSim(iter(id),id)), ' \n'])
            tic
            ep.restart()
            timeRestart(iter(id),id) = toc;
            fprintf(['Restart Time Elapsed: ', num2str(timeRestart(iter(id),id)), ' \n'])
            disp(['Successful reconnection for worker ',num2str(id),'.'])
            if TAKEMEANTIME == 1
                timeSimAvg = timeSim; timeSimAvg(~any(timeSimAvg,2),:) = []; % Empty out zeros for now
                timeSimAvg = mean(timeSimAvg,1); 				 % Average of adjusted array
                timeRestartAvg = timeRestart; timeRestartAvg(~any(timeRestartAvg,2),:) = []; 
                timeRestartAvg = mean(timeRestartAvg,1);
                fprintf(['timeSim Average: ', num2str(timeSimAvg(id)), 's\n']);
                    fprintf(['timeRestart Average: ', num2str(timeRestartAvg(id)), '\n']);
            end 
            tic
        end
    end 
    dayTime = mod(eptime, 86400);  % Current time this simulation day.
    thisHour = floor(dayTime/3600);
    %nextHour = ceil(dayTime/3600);
    %if nextHour == 0, nextHour = 1; end
 
% Mode 0
%------------
% - Off mode used for evaluating NetPenalty Function. (remove)
   if MODE == 0
       
       if coolDownCount == 0 && BUFFERDAY == BUFFERLIMIT % Off for candidate day
           ReheatStatus = zeros(1,5); % Turn off EVERYTHING for the test day.
           CoolingCoil_ON = 0; 
       else  %Tamb = outputs(1); 
        Tzones = outputs(2:7);
        
        %-----------------------
        ReheatStatus = 1.*ones(1,5);   % All reheat on
        CoolingCoil_ON = 1;            % Cooling coil starts on
        deadTol = 0.556;               % New deadzone parameter. Reheat ends and cooling starts 1F below the setpoint.
        for n = 0:1:4
            setPoint = Tset(thisHour+1+24*n);
                 if setPoint > (Tzones(n+1)+deadTol)  % No heating coil service if Tset is high enough
                     ReheatStatus(n+1) = 0; 
                 end
        end
        occupiedHrs = 8:1:19;
        if any(ReheatStatus) == 0 && any(occupiedHrs == thisHour) == 0 % If all heating coils are off, assume no zone should be served at this time
            CoolingCoil_ON = 0;
        end 


            
        
       end 
% Mode 1
%------------
% - Tset behaves as an actual list of 25 set points 
% - Adjust on/off based on zone temperatures

    elseif MODE == 1
        %Tamb = outputs(1); 
        Tzones = outputs(2:7);
        
        %-----------------------
        ReheatStatus = 1.*ones(1,5);   % All reheat on
        CoolingCoil_ON = 1;            % Cooling coil starts on
        deadTol = 0.556;               % New deadzone parameter. Reheat ends and cooling starts 1F below the setpoint.
        for n = 0:1:4
            setPoint = Tset(thisHour+1+24*n);
                 if setPoint > (Tzones(n+1)+deadTol)  % No heating coil service if Tset is high enough
                     ReheatStatus(n+1) = 0; 
                 end
        end
        
        occupiedHrs = 8:1:18;
        if any(ReheatStatus) == 0 && any(occupiedHrs == thisHour) == 0   % If all heating coils are off, assume no zone should be served at this time
            CoolingCoil_ON = 0;
        end 

%-----------------------
   end 
   
  %  SP = [setPoint setPoint HVAC_ON];   % Vector to write.
     SP = [Tset(thisHour+1),Tset(thisHour+1+24*1),Tset(thisHour+1+24*2),Tset(thisHour+1+24*3),Tset(thisHour+1+24*4), ...
           ReheatStatus,CoolingCoil_ON];
       
    % Write out to EnergyPlus.
    ep.write(mlepEncodeRealData(VERNUMBER, 0, (kStep-1)*DELTATIME, SP)); % (kStep-1)*(COOLDOWNCOUNT)?   

    % Save to logdata.
    logData(kStep, :) = horzcat(SP, outputs);   % Check NUMCOLUMNS
    
    if kStep == NUMSTEPS*24*dayFlag        
       if BUFFERDAY == BUFFERLIMIT
            if ~isempty(dayLog)         % Fixes labeling in the structure by counting existing days done.  
                labelOffset = (length(fieldnames(dayLog))+1);
            else
                labelOffset = 1;
            end 
       
            %thisDay = ['Day',num2str(dayFlag-BUFFERLIMIT*labelOffset)];
            thisDay = ['Day',num2str(BUFFERDAY),num2str(coolDownCount)];
            dayLog.(thisDay) = logData(prevkStep:kStep, :); 
            dayFlag = dayFlag+1;
            prevkStep = kStep+1;    % Always mark end of each day. 
            kStep = kStep+1;
            fprintf('Returning.\n')
            return   % Get out to start processing.
       else
             thisDay = ['Day',num2str(BUFFERDAY),num2str(coolDownCount)];
             dayLog.(thisDay) = logData(prevkStep:kStep, :); 
             fprintf('BufferDay = %d / %d\n',BUFFERDAY,BUFFERLIMIT) 
             prevkStep = kStep+1;    % Always mark end of each day. 
             dayFlag = dayFlag+1;
             BUFFERDAY = BUFFERDAY+1;
       end 
    end    
    
    kStep = kStep+1;
  
end

mlepflag
ep.stop;    % STOP THE SIMULATION
thisDay = 'DONE';
% Remove unused entries in logdata
if kStep-1 < MAXSTEPS
    logData((kStep):end,:) = [];
end

end 

% ==========FLAGS==============
% Flag	Description
% +1	Simulation reached end time. (Asking for excess days in this file.)
% 0	    Normal operation.
% -1	Simulation terminated due to an unspecified error.
% -10	Simulation terminated due to an error during the initialization.
% -20	Simulation terminated due to an error during the time integration.

