function [PENNetTotal,PEN_zone,populationCoeff,PelecCost,PenCost] = NetPenalty(Pelec,Tzone,tbegin,tend,interval,div,disableTempPenalty,TambVec,CALIB,cooldownCount,c1,IDFweekday)
global hideFigures
% ================================================
% Constants Box
% ================================================
% Procedure to choosing a new arbitrary WTP curve based on the example:
% 1) Make Tset = whatever setpoints we are currently employing. (Set all to
%    Tideal is typical).
% 2) Disable occupant penalties. (Already negligible if Tset=Tideal)
% 3) MAKE SURE THE COOLING IS ACTUALLY ON
% 4) Simulate and read NetCost. This should be purely the electricity cost
% 5) Therefore, in absence of all cooling, we should be incurring a
%    temperature cost penality (integrated population penalties) of the
%    same amount. Trial and error guess Co until you get approximately the
%    aforementioned NetCost. 

% Modes
%-------------------------------------------------
% If brickWallPenalty = 0, the parabolic temperature penalty is used.
% Otherwise, all temperatures within the range below are treated equally,
% and we are intolerable outside this range. 
% Temperature penalty must be enabled (disableTemperaturePenalty = 0).
    brickWallPenalty = 0;       % "Brick wall mode" with temperature band. 
    Tlow = 66;  Thigh = 76;     % Defines the acceptable tolerance zone for the above.
    
% Plotting Options
%-------------------------------------------------
    doPenaltyPlot = 0;      doZoneTempPlot = 0; 
    doPopulationPlot = 1;   doPENPlot = 1;
    doWTPPlot = 0;
    
% WTP / Penalty Options
%-------------------------------------------------
    %f = 1.5;                     % Arbitrary correctional/adjustment factor (original: 1.5)
    T = [60:.01:100];             % Temperature (F) to plot over.
    %S = 0.0875;                  % Dollars/hour/person for current AC bill. UH ONLY
    %x = 2; Co = x*S;             % Willing to pay x times the current heating bill for Tideal
    
    %Co = 1;           % Tuned 5 Zone: Tideal = 70F during occupied. Tamb = 76.7224.
                         % For the transition penalty, Ct = 5/size(Tzone,1), about 5/287.
                         
    Tideal = 72.5;         % Degrees F, where the occupants pay no temperature penalty. 
    %Tamb = 83;          % Average ambient temperature. Tamb = 83*ones(div*24,1);
    Tamb = mean(TambVec);
    
    if CALIB == 1
        c1 = 1;
    end
    
    C1 = (c1).*ones(length(TambVec),1);
    
    %C1avg = Co*ones(length(Tamb),1)./(Tamb-Tideal).^2;  % C1 = kron(C1,ones(1,size(Tzone,2)));
    %C1 = Co.*ones(length(TambVec),1)./(TambVec-Tideal).^2
  
    
    C1 = kron(C1,ones(1,size(Tzone,2)));

% Population
%-------------------------------------------------

%**************** 5 Zone ****************
% For 4/22.

thisDay = IDFweekday;
dayNames = {'Monday','Tuesday','Wednesday','Thursday','Friday','Saturday','Sunday'};
dayInd = find(strncmpi(thisDay,dayNames,3) == 1, 1,'first');
dayIndSat = 6-dayInd; 

    M = 11; % 1; % Scaling factor. Use 11 for the large office.
    numZones = 5;
    %people = [11 5 11 5 20].*M;            % People located in each zone, based on EP data (small office)
    people = [20,12,20,12,145];             % Large office.
    scheduleOccupy = [0, 8, 11, 12 ,13 ,14 ,18, 19, 24;   % Top row: Until: xx:00 (military time)
                      0, 0, 1,  .8, .4, .8, 1,  .5, 0];   % Bottom row: coefficient (%) of normal population present.
    if any(cooldownCount == [dayIndSat,dayIndSat+1,dayIndSat+7,dayIndSat+1+7]) % Weekend scheduling.
     scheduleOccupy = [0, 8, 9,   11  ,13 ,16 ,18, 19, 24;   % Top row: Until: xx:00 (military time)
                       0, 0, 0.2, 0.3, 0, .3, .2,  .1, 0];   % Bottom row: coefficient (%) of normal population present.
    end 
                  
                  
                  
%**************** 60 Zone ****************
% For 8/22
% Note we really export 9 unique lumped zones.
% 1440 m^2/60 zones = 24 m^2/zone
% 24 m^2/zone * 1/9.29 persons/m^2  ~ 2.58 persons/zone

%     numZones = 9;
%     people = 24/9.29*ones(1,9); 
%     scheduleOccupy = [0, 4, 5, 6, 7, 11, 13, 17, 18, 19, 20, 23, 24;
%                       0, 0, 0.05, 0.1, 0.25, 0.65, 0.6, 0.65, 0.4, 0.25, 0.1, 0.05, 0];
% ================================================

    if hideFigures == 1
        doPenaltyPlot = 0;      doZoneTempPlot = 0; 
        doPopulationPlot = 0;   doPENPlot = 0;
        doWTPPlot = 0;
    end 
    
%% Population Scheduling
% Using single population schedule for all zones.            
    occupyCoeff1 = genOccupyCoeff(scheduleOccupy,div);    % Generate a schedule's coefficients using this function.
    occupyCoeffZones = repmat(occupyCoeff1,[1,numZones]);        % Population coefficients for all zones (duplicated 1 zone)
    populationCoeff = repmat(people,24*div,1) .*occupyCoeffZones;   % Multiply element-wise into actual population values
    if doPopulationPlot == 1
        figure(2), plotDate(populationCoeff,tbegin,tend,[1,1,1,60],'-*'); ylabel('Population'), xlabel('Time (5 Min Intervals)'), grid on, title('Population over 24hrs')
    end 
    
%% Willingness to Pay

% WTP/PEN: Raw concave parabolic WTP/PEN curve.
%   - Comfort relationship is assumed to be quadratic
%   - Quantities are based on the entire population willing to pay up to
%     twice the normal price for comfort. 
%   - See Constants Box for more info
% WTP/PEN_zone: Functions evaluated based on individual Tzone temperatures,
%               resulting a curve for each zone.
% PENZoneTotal: Combined zone penalty curve.
% PelecInterp:  Interpolated HVAC price curve (based on HECO prices) since
%               prices are supplied hourly. Used to properly add price vectors.
% PENNetTotal: PENZoneTotal + PelecInterp (master penalty curve)

   % WTP = WillToPay(T,Co,C1avg,Tideal);    WTP = WTP*sum(people);
   % PEN  = PenaltyPay(T,C1,Tideal);     PEN =  PEN*sum(people);
    
    if doWTPPlot == 1
        figure(3)
        plot(T,WTP,'b'), grid on, title('"Willingness to Pay"'), hold on
        plot(T,PEN,'r'), grid on, title('"Penalty"'), ylim([0,max(WTP)]), xlabel('Desired Temperature (F)'), ylabel('$/hour/day')
    end

    % Should have 288 points: 5 min sampling, 12 per hour, 24 hours. 
    %WTP_zone = WillToPay(Tzone,Co,C1,Tideal).* populationCoeff;%-0.01*abs(repmat(TsetList,1,size(Tzone,2))-Tideal*ones(length(TsetList),size(Tzone,2)));        
    PEN_zone = PenaltyPay(Tzone,C1,Tideal).* populationCoeff;%+0.01*abs(repmat(TsetList,1,size(Tzone,2))-Tideal*ones(length(TsetList),size(Tzone,2)));        %+0.1*abs(TsetList-Tideal);
    
    % This is for the "brick wall" mode.
    if brickWallPenalty == 1
        PEN_zone = 0;                               % All comfort zone temperatures are treated as equally tolerable.
        TzoneComp = Tzone.*(populationCoeff > 0);   % Find populated hours only...
        [nzRows,~] = find(TzoneComp > 0);           % Return only those Tzone rows
        nzRows = unique(nzRows);                    % No duplicates
        TzoneComp = TzoneComp(nzRows,:);
        
        % Reshape the populated zone temperatures into hours for Tavg per
        % hour for all zones. 12 is the timesteps in each hour (5 min).
        finalZoneComp = [];
        for i = 1:size(TzoneComp,2)
           TsingleZone = TzoneComp(:,i);
           TsingleZone = reshape(TsingleZone,12,size(TsingleZone,1)/12);
           TsingleZone = [mean(TsingleZone(3:end,:),1)]; % Take over each hour
           finalZoneComp = [finalZoneComp,reshape(TsingleZone,numel(TsingleZone),1)];
        end
        TzoneComp = finalZoneComp;
        if any(any(TzoneComp > Thigh | TzoneComp < Tlow));  % If any populated zone temperatures are outside the margins,
            PEN_zone = (10e3)*ones(288,5); % Arbitrary brick wall term is added to the energy cost. This is technically in Wh.
            fprintf(2,'(Hit Brick Wall)')
        else
            PEN_zone = zeros(288,5);
        end
    end 
    
    if doPenaltyPlot == 1
        figure(4), plotZones(PEN_zone,tbegin,tend,interval,numZones,'Comfort Penalty ($/hr/day)', doPenaltyPlot)
    end 
    if doZoneTempPlot == 1
        figure(5), plotZones(Tzone,tbegin,tend,interval,numZones,'Temperature (F)', doZoneTempPlot)
    end 
    
    PENZoneTotal = sum(PEN_zone,2);  %Pelec = [Pelec,Pelec(1)];  
    
    if disableTempPenalty == 1
        PENZoneTotal = PENZoneTotal == -1; % Disable temp penalties by setting them all to zero.
        fprintf('(PENALTY DISABLED)')
    end 
   % PelecInterp = interp1(1:25,Pelec,1:24/(length(PENZoneTotal)-1):25);
   
   %[l,w] = size(Pelec)
    PelecInterp = Pelec;  % Must change marked place in plotHVAC
    %PelecInterp = kron(Pelec,ones(1,div))';
    
    PelecCost = trapz(linspace(1/div,24,24*div),PelecInterp,1);  % Electricity cos
    PenCost = trapz(linspace(1/div,24,24*div),PENZoneTotal,1);   % Penalty cost 
    PENNetTotal = PelecInterp+PENZoneTotal;                      % Summation of both costs
    
    
    % Plotting functions
    if doPENPlot == 1
        figure(4)
        subplot(3,1,1)
        plotDate(PelecInterp,tbegin,tend,interval,'-b.'); ylabel('$/hr Penalty'), xlabel('Time (30 Min Intervals)'), grid on, title('HVAC Price Interpolated'),hold all
        rotateticklabel(gca);
        subplot(3,1,2)
        plotDate(PENZoneTotal,tbegin,tend,interval,'-m.'); ylabel('$/hr Penalty'), xlabel('Time (30 Min Intervals)'), grid on, title('Temperature Penalty')
        rotateticklabel(gca);
        subplot(3,1,3)
        plotDate(PENNetTotal,tbegin,tend,interval,'-r.'); ylabel('$/hr Penalty'), xlabel('Time (30 Min Intervals)'), grid on, title('Summed Penalties')
        rotateticklabel(gca);
    end 

end 


%% Population Coefficients
% Generates population coefficients to multiply into an occupancy matrix.
% Expects something like the following as its input:

% scheduleOccupy = [0, 8, 11, 12 ,13 ,14 ,18, 19, 24; 
%                   0, 0, 1,  .8, .4, .8, 1,  .5, 0];
% Top Row: "Until:" parameters as listed in EP. Needs to have 0 and 24 at
%           the beginning and end for this to work properly
% Bottom Row: Coefficients (0-100%) occupancy to model population movement.

function occupyCoeff = genOccupyCoeff(scheduleOccupy,div)
    FIVEMINOFFSET = 1;  % We start at 00:05:00 in EP (5 min after midnight), so there should be one LESS coefficient for this interval.
    occupyCoeff = zeros(24*div,1);          % Should have 24 hrs * 'n' points/hr coefficients.
    divCount = 1;                           % Marker to begin writing in new coefficients into master vector
    for i = 2:size(scheduleOccupy,2)        % Use 2 to avoid referencing index 0
        untilDiv = ones((scheduleOccupy(1,i)-scheduleOccupy(1,i-1))*div-FIVEMINOFFSET,1); % Hourly Time Displacement * points/hr 
        occupyCoeff(divCount:divCount+length(untilDiv)-1) = kron(scheduleOccupy(2,i),untilDiv);
        % Use kron to replicate coefficient length(untilDiv) times, then
        % merge with previous coefficients using the assignment marker
        divCount = divCount + length(untilDiv);  % Marker should be updated AFTER calculation. 
        FIVEMINOFFSET = 0;  % After the very first interval (00:00:00), resume having 12 divs per hour.
    end 
        %hourValues = 1/12:1/12:24;     % For debug testing, make sure times line up.
        %occupyCoeff = [hourValues',occupyCoeff];
end 


%% WTP / PP Base Functions

function [val] = WillToPay(T,Co,C1,Tideal)
        val = -C1.*(T-Tideal*ones(size(T,1),size(T,2))).^2+Co;  
end 

function [val] = PenaltyPay(T,C1,Tideal)
   % if size(T,1) == 1;
        val = C1.*(T-Tideal*ones(size(T,1),size(T,2))).^2;
   % else
   %     val = C1.*(T-Tideal).^2 + C1.*[zeros(1,5); T(2:end,:)-T(1:end-1,:) ].^2; % Use 1/25 to punish a 5F change in 5 min (1 timestep) with 1$. 
    %end 
end 

%% END



