function [NetCostTOTAL,PelecTOTAL,PdiscTOTAL,WattTOTAL] = mainEP(Tset,MODE,CALIB,doy,numMultiDays,testMODE2,c1Str,searchname,dateObj,IDFNAME)
    % These globals are to improve efficiency of iterations and to return
    % extra information done throughout Nelder Mead as we converge. 
    global hideFigures  % Do not affect computations
    global iter, global cost, global ep
    global calibArr
    global outputTable
    
    % Options
    %---------------------
    customLambda = 0;    % Override with custom system lambda (overwrites cost(DOY)) (Obsolete)
    numWorkers = 15;     % How many workers are there? Ignored if MODE = 0.
    COOLDOWNDAYS = 13;    % 3 or 13; % How many cooldown days? (0 means only test day).
    stopFigures = 1;     % For single run mode: hide figures?
    incrementDOY = 1;    % Should non-custom lambda be incremented? 
    doPause = 0;         % For single runs only; should be 0 for multirun
    IDFNAME = [IDFNAME,'.idf'];
    WEATHERNAME = 'USA_IL_Chicago-OHare.Intl.AP.725300_TMY3.epw';
    %WEATHERNAME = 'USA_HI_Honolulu.Intl.AP.911820_TMY3.epw';  % Doesn't open properly? Make sure to copy into each directory
    saveTestDirectoryName = ['HorizonLogs_',c1Str];
    
                                % CHANGE REPEATS IN THE IDF to 3
    pickUpResults = 0;          % If true, adds 1 to COOLDOWNDAYS to form an extra day to terminate on. (Obsolete)
    disableTempPenalty = 0;     % Disable temperature penalty entirely. 
    disableTransPenalty = 1;    % Disable transition temperature penalty entirely.  (Obsolete)
    stepPref = 0;               % Do step temperature preference (overwrites Tset). (Obsolete)
    showTset = 0;               % Show setpoints in figure 5? (Overrides CALIB and stepPref). (Obsolete)
    useTerminationHorizon = 1;  % Include cool down days in total cost too?
    useSteadyState = 0;         % Use the latest summed NetCost value? (use with useTerminationHorizon = 0).
                                % Also remember to use non-incrementing days in the IDF. (Obsolete)
    if strcmp(c1Str,'10e-3')
        disableTempPenalty = 1;
        useTerminationHorizon = 1;
        saveTestDirectoryName = ['HorizonLogs_Corbin_WithTerm_',c1Str];
    elseif strcmp(c1Str,'20e-3')
        disableTempPenalty = 1;
        useTerminationHorizon = 0;
        saveTestDirectoryName = ['HorizonLogs_Corbin_NoTerm',c1Str];
    end 
    mkdir(['/home/atizawa/lus/Window_Shifting_Logs/',saveTestDirectoryName])                          
                                
    c1 = str2double(c1Str); 
    TsetHistory = readTsetHistory(searchname,numMultiDays); % 2) Assign historical policy
    
    outputTable = zeros(COOLDOWNDAYS+1,7);
    [directoryName,tableName,~] = fileparts(searchname);
    [~,IDFweekday] = weekday(dateObj,'long');

    % Temperature preferences/options
    %---------------------
    TbaseF = 75;     % Tbase in F, (not absolute temperature) > Tideal
    dT = (1/1.8)/2;  % 1/2F deviation in C, used for C1 calibration.
    TprefL = 75;     % stepPref method low bound (F)
    TprefH = 79;     % stepPref method high bound (F)
    numPoints = 5;  % Number of temperature (and cost) ranges to make.
    
    %--------------------------------------------------
   
    
    NetCostTOTAL = 0;
    PelecTOTAL = 0; PdiscTOTAL = 0;
    Tbase = convtemp(TbaseF,'F','C');   % Temperature in Celsius that we are WTP (Tideal + n)
    
    % Overrides
    
    if MODE == 1 || MODE == 2 && testMODE2 == 0, stopFigures = 1; end
    if MODE == 0, numWorkers = 1; end 
    if testMODE2 == 1, useTerminationHorizon = 1; disableTempPenalty = 0; numWorkers = 1; end  %stopFigures = 0;
    if useSteadyState == 1, useTerminationHorizon = 0; incrementDOY = 0; end  
    
    id = get(getCurrentTask,'ID');
    
    
        
    % Derivative zone scaling and translation code
    %--------------------------------------------------
    
     if MODE ~= 2
        % For modes 0 and 1, change the 64D vector into a 120 length
        % vector (change the format).
        PSOTset = Tset;
        for z = 0:4
            startInd = 9; endInd = 19;
            duration = endInd-startInd+1;
            lowInd = startInd+duration*z; highInd = endInd+duration*z; 
            Tset(1+24*z:24*(z+1)) = [PSOTset(1:8),PSOTset(lowInd:highInd),PSOTset(end).*ones(1,24-endInd)];
        end
    end 
    
    if MODE == 2
       multiDayTset = Tset; 
    end
    
    useSetPoints = Tset; 
    % From here on, Tset is 24*5 dimensional. 
    if pickUpResults == 1, COOLDOWNDAYS = COOLDOWNDAYS+1; end
    doyInt = doy; % Start with this doy
    for cooldownCount = 0:1:COOLDOWNDAYS
        if incrementDOY == 1
            doy = mod(doyInt+cooldownCount-1,365)+1; % Increment DOY each time, accounting for year-wrapping
        end
        % Tset code for multi day support. Convert from 98 -> 120 length
        % vector, based on occupied/unoccupied period.
        if MODE == 2 && cooldownCount < numMultiDays
           Tset = convForMODE2(cooldownCount,multiDayTset);
        end
    
        
        myDefaultScheduleTset = [];
        
        
        for j = 1:numMultiDays
           myDefaultScheduleTset(j,:) = convForMODE2(0,TsetHistory(j,:));   % Pretend it is the first day         
        end

        % This is the default schedule for the cooldown days. 
        if cooldownCount >= numMultiDays %(cooldownCount > 0  && MODE ~= 2) || (MODE == 2 && cooldownCount >= numMultiDays)
            Tset = myDefaultScheduleTset(cooldownCount-numMultiDays+1,:); % Termination period application
        end 

        if isempty(iter) % It's our first try...
            fprintf('**** First iteration ****\n')
            iter = ones(1,numWorkers);
        end
        
        if isempty(cost)
           cost = importLambda();   % Import lambda for price data
           if customLambda == 1  % Overwrite with custom if allowed.
               % Custom lambda
               fprintf(2,'Custom system lambda being used.\n')
               M = 120;  % Maximum $/MWh (offset removed)
               O = 30;   % $/MWh offset; minimum cost for some daily hour.
               t = -.3:(.2+.3)/23:.2;
               cost = sinc(2*pi*t);
               cost = cost-min(cost); 
               cost = cost./norm(cost,Inf); 
               cost = cost*M+O;
               
               cost = [43.9, 35.26,30.49,27.85,27.67,31.88,37.61,48.78,59.14,72.77,85.48,97.35,106.13,116.69,128.9,139,143.87,131.63,113.09,103.26,101.8,95.47,75.48,58.11];
   
               cost = kron(cost,ones(365,1));  % DOY is irrelevant.
           end 
        end 
        
        dayFlag = 1; kStep = 1; prevkStep = 1; logData = [];  
        
        if cooldownCount == 0
            dayLog = struct(); % Creation of structure for simulation history this iteration.
        end 
        
        if stepPref == 1 && cooldownCount == 0
            if MODE == 2
               error('No support for step temperature preference with multidays yet. Terminating.')
            end
            
            fprintf(2,'Using step temperature policy assignment.\n')
            Tpref = linspace(TprefL,TprefH,numPoints); % Single degree F temperature steps
            todaysPrices = cost(doy,:);                    % For this day's prices,
            costSDim = sort(reshape(cost,numel(cost),1));
            
            percentiles = linspace(0,100,numPoints+1);
            fprintf('Percentile rankings for annual system lambda:\n')
            costBounds = prctile(costSDim,percentiles(2:end-1)); % get the boundaries
            costBoundsL = [0,costBounds];
            costBoundsH = [costBounds,Inf];
            for hr = 1:1:24
               costThisHr = cost(doy,hr);
               ind = find((costBoundsH >= costThisHr) & (costBoundsL < costThisHr));
               Tset(hr) = Tpref(ind);  % Assign the hourly Tset values.
            end
            
            disp(percentiles(2:end-1))
            disp(costBounds)

            fprintf('Assigned the following 24-hr all-zone Tset policy:\n')
            disp(Tset(1:24))
            fprintf('Tset range: '),disp(Tpref),fprintf('\n')
            disp(cost(doy,:))
            Tset = repmat(Tset(1:24),1,5);  % Overwrite existing Tset.
            Tset = convtemp(Tset,'F','C');
            %pause()
        end
        
        %% Options
        %http://blogs.mathworks.com/loren/2009/10/02/using-parfor-loops-getting-up-and-running
        %======================================================================
        % OPTIONS
        %======================================================================
        if MODE == 0 || MODE == 2
            if CALIB == 0
                hideFigures = stopFigures;                 % For more: toggle Figures 4-7 in NetPenalty.m. 
            else
                hideFigures = 1;
            end
        else
            hideFigures = 1;
        end                                  
                                             % Only called here, plotDate and NetPenalty.m
            tbegin = '22-Jun-2013 00:00:00';    % Referenced and manipulated in plotHVAC and plotTemperature.
            tend = '22-Jun-2013 23:55:00';      % Change the date only (not time).
            div = 12;    % Divisons per hour. For five minute intervals, 12 samples per hour.
                         % 60/div should match the timestep value in the IDF. 
           % doy = DOY;   % doy: Corresponds to tbegin and tend day. Refer to Lambda doc: this value is also the row number 
                         % of the correct date minus 1 (killed the header when importing the data). 
                         % cost(n,:) corresponds to the day of the year in 2013
                         % (out of 365). doy = 233 is 8/21. 

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


        %% Simulation and Setup

%             if isempty(id) && MODE ~= 1    % When the worker terminates...
%                 NetCostTOTAL = 99999;      % Something really large such that PSO doesn't think NetCostTOTAL = 0. 
%                 return
%             else
            if isempty(id) && MODE ~= 1  % For usage without the parallel pool.
                id = 1; 
            end 
            if CALIB == 1 &&  cooldownCount == 0
                fprintf(2,'########################################\n')
                fprintf(2,'Doing calibration sequence. [%d/5]\n',iter(id))
                fprintf(2,'########################################\n')              
                            
                Tset = [35*ones(1,24), zeros(1,4), zeros(1,4)];  % 24 setpoints with 8 tweaking points for remaining 4 zones.
                if iter(id) == 1                          % Test 1: Pelec+Pmin (left point)
                    disableTempPenalty = 1; 
                    Tset(9:19) = (Tbase-dT)*ones(1,11);  
                elseif iter(id) == 2                      % Test 2: Pmin only
                    disableTempPenalty = 1;               
                    Tset(9:19) = 32*ones(1,11); 
                elseif iter(id) == 3                      % Test 3: Pelec+Pmin (right point) 
                    disableTempPenalty = 1;
                    Tset(9:19) = (Tbase+dT)*ones(1,11);
                elseif iter(id) == 4                      % Test 4: Pelec+Pmin+Pdisc (left point)
                    disableTempPenalty = 0;
                    Tset(9:19) = (Tbase-dT)*ones(1,11);
                elseif iter(id) == 5                      % Test 5: Pelec+Pmin+Pdisc (right point)
                    disableTempPenalty = 0;
                    Tset(9:19) = (Tbase+dT)*ones(1,11);
                end
                    Tset = repmat(Tset(1:24),1,5);
           end                
            
            
            %   [Tamb, Tzone, PHVAC] = simEP(Tset);       % Run simulation and import

                USECHECKPOINT = 0;             % Use checkpoint?
                TIMESTEP = 60/div;
                BUFFERLIMIT = 3;               % DISABLE IF USING CHECKPOINT BY PUTTING 0. How many days to buffer with default schedule? 
                %NUMCOLUMNS = 3+7;          % Should be number returned parameters. Should be length(SP) + num returned parameters
                NUMCOLUMNS = 18+9;
                NUMDAYS = 9999;           % must match number of times to repeat in IDF
                NUMSTEPS = 60/TIMESTEP;  % Matches numsteps/hr as stated in IDF exactly
                MAXSTEPS = NUMSTEPS*24*NUMDAYS*(BUFFERLIMIT+1);  

                if iter(id) == 1 && cooldownCount == 0 
                    fprintf('Beginning first simulation. Warming up (may take a while).\n')
                    %pause((id-1)*5)
                    ep = mlepProcess;
                    move2EPdirectory(IDFNAME,c1Str);
                    if USECHECKPOINT == 1
                        %ep.program = ['plink -ssh 192.168.56.101 -l Izawa -t -R {port}:localhost:{port} ' ....
                        %          'cd /home/osboxes/bcvtb/examples/5ZoneTsetTest && dmtcp_launch --new-coordinator --no-gzip python ./ep_worker {input_file} {weather_file} {port}"'];

                        portFile = ['portWorker', num2str(id)];                                % File name for this worker's DMTCP port.
                        portCmd = ['dmtcp_coordinator -p 0 --ckptdir /home/atizawa/5ZoneTsetTest --port-file ',portFile, ' &'];   % Start coordinator and note the port.
                        system(portCmd,'-echo');
                        fprintf(['Started coordinator with:', portCmd,'\n'])
                        % ERROR BELOW: Generate valid file identifier. Insert
                        % pause statement before opening?

                        retryLimit = 5;
                        for tryOpen = 1:1:(retryLimit-1)
                            pause(2); % Wait a while; the coordinator needs time to write the file.
                            fid = fopen(portFile,'r'); % Try to get the port number for this worker's coordinator.
                            if fid == -1
                                fprintf(2,['WORKER ', num2str(id), ' failed to open ',portfile,'. Retrying in 1 second. (ATTEMPT #%d',tryOpen,')\n'])  
                            else 
                                dmtcp_port = fscanf(fid,'%d'); fclose(fid); % Get the DMTCP port to feed to Python.
                                fprintf(['WORKER ', num2str(id),': retreived DMTCP port ', num2str(dmtcp_port),'\n'])
                                break
                            end
                        end 
                        if tryOpen == retryLimit, fprintf(2,'Failed to open after %d attempts. breaking.\n',retryLimit),end

                        ep.program = ['ssh atizawa@redr.eng.hawaii.edu -T -R {port}:localhost:{port} ' ...
                                      '"cd /home/atizawa/5ZoneTsetTest && dmtcp_launch -p ', num2str(dmtcp_port), ' --no-gzip python ./ep_worker {input_file} {weather_file} {port} ', num2str(dmtcp_port), '"'];
                        % The dmtcp port number is passed in twice since we'll
                        % need it to reconnect later.
                        ep.acceptTimeout = 200000; % Change this if acceptTimeOut appears. Was  originally  10000
                        else
                           %ep.program = ['plink -ssh redr.eng.hawaii.edu -l atizawa -pw toshio3566 -t -R {port}:localhost:{port} ' ...
                           %              '"cd /home/atizawa/5ZoneTsetTest && export BCVTB_HOME=/home/atizawa/5ZoneTsetTest/bcvtb && runenergyplus SmallOffice2004TsetTest.idf USA_HI_Honolulu.911820_TMY2.epw"']; % No Python helper server
                           ep.program = ['export BCVTB_HOME=/home/atizawa/apps/RefLargeOfficeBuilding/bcvtb && /home/atizawa/apps/EnergyPlus-8-3-0/runenergyplus']; % No Python helper server

                            ep.acceptTimeout = 200000; % Change this if acceptTimeOut appears. Was  originally  10000
                        end 

                    ep.arguments = {IDFNAME, WEATHERNAME};
                    fprintf('Waiting for EnergyPlus socket...\n')
                    [status, msg] = ep.start;  
                    if status ~= 0, error('Could not start EnergyPlus: %s.', msg); end
                    fprintf('Got it!\n')
                    logData = zeros(MAXSTEPS, NUMCOLUMNS);      % Holds all output data.
                    cd '/home/atizawa/apps/Cost-Analysis-Model' % USE IF ON LOCALHOST
                end 
            
                [dayLog,thisDay,kStep] = doMLEPSim(Tset,dayLog,ep,60/div,dayFlag,kStep,prevkStep,logData,iter,cooldownCount,COOLDOWNDAYS,Tbase,numMultiDays,myDefaultScheduleTset);
                fprintf('MLEPSim done.\n')
                if strcmp(thisDay,'DONE')
                   if pickUpResults == 0 
                       fprintf(2,'Done simulating due to EP day limit at "%s".\n',thisDay)
                       fprintf(2,'Check settings in EP or in doMLEP.\n')
                       fprintf(2,'If 6 days were set in the IDF runperiod, then csv outputs were provided in EP_Worker1 with proper termination.')
                       error('Check if results converged, or add more iterations. Check variable.cfg. Prematurely ended fminsearch.')
                   else
                       fprintf(2,'Fetching results and terminating NOW.\n')
                       fprintf(2,[IDFNAME,' has ESO results available in EP_Worker_1/Output.\n'])
                       error('Breaking. If this was not intentional, set pickUpResults = 0.')
                   end 
                end 
%                 Tamb = dayLog.(thisDay)(:,4);
%                 Tzone = dayLog.(thisDay)(:,5:9);%end-1);
%                 PHVAC = dayLog.(thisDay)(:,10);
%                 TsetList = dayLog.(thisDay)(:,1);   % For plotting purposes only
%                 TsetList = convtemp(TsetList,'C','F');  Tamb = convtemp(Tamb,'C','F'); 
%                 Tzone = convtemp(Tzone,'C','F');
%                 CoolingPumpAvailSched = dayLog.(thisDay)(:,3);

                Tamb = dayLog.(thisDay)(:,11+1);
                Tzone = dayLog.(thisDay)(:,11+2:11+6);%end-1);
                PHVAC = dayLog.(thisDay)(:,11+7);
                TsetList = dayLog.(thisDay)(:,1);   % For plotting purposes only
                TsetList = convtemp(TsetList,'C','F');  Tamb = convtemp(Tamb,'C','F'); 
                Tzone = convtemp(Tzone,'C','F');
                CoolingPumpAvailSched = dayLog.(thisDay)(:,11);

               % Plot temperature setpoints (ignore cost, comment out)
                 if cooldownCount == 0 && showTset == 1
                     useSetPoints = reshape(useSetPoints,24,5);
                     Tzone = kron(useSetPoints,ones(12,1)); 
                     Tzone = convtemp(Tzone,'C','F');
                 end 
                
                
        %% Computations

            if hideFigures == 1
                close all, set(0,'DefaultFigureVisible','off')
            else
                set(0,'DefaultFigureVisible','on')
            end 
            if cooldownCount == 0
                fprintf(['\n============================= Iteration ', num2str(iter(id)), ' | Worker ', num2str(id),'. =============================\n'])
                fprintf('Running computations for:\n')
                fprintf([tbegin,'\n',tend,'\n'])
            end 
            tstep = [1 1 1 30];
            [Pelec, wattHoursConsumed] = plotHVAC(PHVAC,cost,doy,tstep,div,hideFigures);          fprintf('...')            
            [PENNetTotal,PEN_zone,populationCoeff, PelecCost,PenCost] = NetPenalty(Pelec,Tzone,tbegin,tend,tstep,div,disableTempPenalty,Tamb,CALIB,cooldownCount,c1,IDFweekday);      fprintf('...')
            NetCost = trapz(linspace(1/div,24,24*div),PENNetTotal,1);                             fprintf('...')
            plotTemperature(Tamb,TsetList,Tzone,tbegin,tend,hideFigures,CoolingPumpAvailSched);   fprintf('...\n')
            
            transitionPenalty = getTransitionPenalty(Tzone,populationCoeff);
            if disableTransPenalty == 1, transitionPenalty = 0; end 
            NetCost = NetCost + transitionPenalty;
        %% Output
        if cooldownCount == 0, Tset; fprintf('Writing out values.\n'), end 
        fprintf('For DOY = %d:\n',doy)
        fprintf('Electricity Cost of Day (%d): $%f\n',cooldownCount, PelecCost)
        fprintf('Discomfort Cost of Day (%d): $%f\n',cooldownCount, PenCost)     
        fprintf('Daily Cost of Day (%d): $%f\n',cooldownCount, NetCost)
        fprintf('Transition Penalty (%d): $%f\n',cooldownCount,transitionPenalty);
        wattHoursConsumed = kron(wattHoursConsumed,ones(1,div))';  % Reported back on per hour basis, go to timestep
        wattHoursConsumed = trapz(linspace(1/div,24,24*div),wattHoursConsumed,1);                          
        fprintf('WattHours Consumed on (%d): %f Wh\n',cooldownCount,wattHoursConsumed)
        
       checkCond(1) = (cooldownCount == 0);
       checkCond(2) = (useTerminationHorizon == 1);
       checkCond(3) = (cooldownCount < numMultiDays);
       
        if any(checkCond)
                PelecTOTAL = PelecTOTAL + PelecCost;
                PdiscTOTAL = PdiscTOTAL + PenCost; 
                NetCostTOTAL = NetCostTOTAL + NetCost;
                
                if useSteadyState == 1
                    NetCostTOTAL = NetCost;
                end 
        else
            fprintf(2,'Ignoring termination horizon %d out of %d. (numMultiDays = %d)\n',cooldownCount,COOLDOWNDAYS,numMultiDays)
        end
            fprintf('TOTAL Pelec: $%f\n',PelecTOTAL)
            fprintf('TOTAL Pdisc: $%f\n',PdiscTOTAL)
            fprintf('TOTAL NetCost: $%f\n',NetCostTOTAL)
            if MODE == 0 || testMODE2 == 1
                outputTable(cooldownCount+1,:) = [PelecCost,PenCost,NetCost,wattHoursConsumed,PelecTOTAL,PdiscTOTAL,NetCostTOTAL];
            end 
            
        if hideFigures == 0   % Plot saving code.
            figureMainDir = [pwd,'/Plots/'];
            fullName =  [figureMainDir,saveTestDirectoryName,'_',tableName,'/Planned_',datestr(dateObj),'_DOY',num2str(doyInt)];
            if exist(fullName,'dir') ~= 2, mkdir(fullName); disp(['Made directory: ' fullName]); end 
            for k = 1:5
                baseFileName = sprintf('fig%d_cooldown%d_DOY%d',k,cooldownCount,doy);
                saveas(figure(k),[fullName,'/',baseFileName,'.fig']);
            end 
            fprintf(['Saved plots to ',pwd,'/Plots/\n'])
            if cooldownCount ~= COOLDOWNDAYS, close all; end % Show last one only.
        end 
    end  % End cooldownCount.
    
    if CALIB == 1
        i = double(iter(id));
        calibArr(i) = NetCostTOTAL;
        if i == 5  % Time to calculate C1 value. 
            PelecL = calibArr(1) - calibArr(2);
            PelecR = calibArr(3) - calibArr(2); 
            PdiscL = calibArr(4) - calibArr(1);
            PdiscR = calibArr(5) - calibArr(3);
            calibC1 = (PelecL-PelecR)/(PdiscR-PdiscL);
        end
    end
    
    
    iter(id) = iter(id)+1; 
    fprintf('Finished iteration.\n')
    
    % Write to workspace for later viewing for single trial.
    if MODE == 0 || (MODE == 2 && testMODE2 == 1)
        dayLog
        fullName = [directoryName,'/',saveTestDirectoryName];
        if exist(fullName,'dir') ~= 2, mkdir(fullName); disp(['Made directory: ' fullName]); end 
        assignin('base','dayLog',dayLog);
        tableName = [tableName,'_doyInt',num2str(doyInt),'.xlsx'];
        fprintf('Saved dayLog. Exporting %s.\n',tableName)
        xlswrite([fullName,'/',tableName],outputTable)
        fprintf('Exporting dayLog')
        writeDayLogStruct2CSV(dayLog,fullName,dateObj,doyInt)
        
        if CALIB == 0
            if doPause == 1
                plot5DayTempScript(dayLog)
                fprintf('Press to continue.\n')
                pause()
            else
                PelecTOTAL = outputTable(1,1);
                PdiscTOTAL = outputTable(1,2);
                NetCostTOTAL = outputTable(1,3);
                WattTOTAL = outputTable(1,4);
                return
            end
        elseif CALIB == 1 && iter(id) == 6 % Post increment iter(id)
            fprintf('**************************************\n')
            fprintf('Finished calibration.\n')
            fprintf('TbaseF = %d F; dT = %f delta C (%f delta F)\n',TbaseF,dT,dT*1.8);
            fprintf('Resulting C1 in ($/person/(delta C)^2/timestep):\n')
            format short eng
            disp(calibC1)
            fprintf('Calibration array:\n')
            disp(calibArr)
            format short
            fprintf(2,'Done with calibration.\n')
            return
        elseif CALIB == 1 && iter(id) > 6
            error('Calibration error; did not terminate properly. Returning.')
        end
    end 
end 

function transitionPenalty = getTransitionPenalty(Tzone,populationCoeff)

    populationCoeffMask = populationCoeff > 0; % Only get Tset for where people are present.
    
    Tzone = Tzone .* populationCoeffMask;
    Tzone = convtemp(Tzone,'F','C');      % Tzone was passed in as F, not C.
    Ct = 5/size(Tzone,1);                 % Coefficient for penalty in $/delta C^2/person. 
    Tdiff = abs(diff(Tzone,1)).^1;        % Should be length 288-1 = 287
    TdiffZonePenalty = (Ct.*ones(size(Tdiff,1),size(Tdiff,2))).*Tdiff.*populationCoeff(1:end-1,:);
    % Apply the coefficient, then scale penalties according to how the zones are populated.
    
    transitionPenalty = sum(sum(TdiffZonePenalty)); 

end 
