% *****      Example 3.4 - ON / OFF Temperature Control     *****
% *****                  from Marlin 1994                   *****
% COPYRIGHT McMASTER UNIVERSITY, 1994
% WRITTEN BY : DAVE JEWELL
% FILENAME : EXAM3_4.M
% LAST UPDATED : NOVEMBER 1994
% VERSION : 1.0
% DESCRIPTION OF FILE :
% THIS PROGRAM SIMULATES AN ON/OFF TEMPERATURE CONTROLLER
% AS DISCUSSED IN THE TEXT AND CARRIES OUT QUESTION 3-9
% ***************************************************************

clear;                       % CLEARS WORKSPACE
clf;                         % CLEARS GRAPH
axis('auto');                % RESETS THE GRAPH AXIS
clc;                         % CLEARS SCREEN

% ***************************************************************
% PROCESS DATA
% ***************************************************************
Cv = 0.171;                  % HEAT CAPACITY OF AIR (cal/g/degC)
rho = 1190.0;                % DENSITY OF AIR (g/m^3)
V = 75.0;                    % ROOM VOLUME (m^3)
UA = 45.0E3;                 % HEAT TRANSFER COEFFICIENT (cal/degC/hour)
Qfon = 1.5E6;                % HEAT ENERGY APPLIED BY FURNACE (cal/hour)
tauQ = 0.0;                  % TIME CONSTANT BETWEEN FURNACE AND ROOM (hour)
Tupper = 23.0;               % UPPER TEMPERATURE CONSTRAINT (degC)
Tlower = 17.0;               % LOWER TEMPERATURE CONSTRAINT (degC)
Ta = 10.0;                   % AMBIENT TEMPERATURE (degC)
Tinit = 20.0;                % INTIAL TEMPERATURE (degC)
Cpwood = 0.5;                % HEAT CAPACITY OF FURNITURE IN ROOM (cal/g/degC)
Mwood = 000.0E3;             % MASS OF FURNITURE (g)
Fair = 0.0;                  % VOLUMETRIC FLOWRATE OF AIR (m^3/hour)

% ***************************************************************
% SIMULATION DATA
% ***************************************************************
tstart = 0.0;                % TIME TO BEGIN SIMULATION (hours)
tend = 2.0;                  % TIME THE SIMULATION ENDS (hours)
delt = 0.002;                % EULER INTEGRATION STEP SIZE (hours)
t = tstart:delt:tend;        % TIME VECTOR USED FOR PLOTTING
n = round((tend - tstart)/delt)+1; % NUMBER OF INTEGRATION STEPS

% ***************************************************************
% INITIAL CONDITIONS AND STORAGE VECTORS
% ***************************************************************
T = zeros(1,n);
Qh = zeros(1,n);             % INITIALIZATION OF STORAGE VECTORS
Qf = zeros(1,n);


T(1) = Tinit;                % INITIALIZATION OF ROOM TEMPERATURE
Qh(1) = 0.0;                 % INITIALIZATION OF HEAT IN ROOM
Qf(1) = 0.0;                 % INITIALIZATION OF FURNACE INPUT

Tmax = ones(1,n).*Tupper;    % VECTOR OF UPPER TEMPERATURE LIMIT
Tmin = ones(1,n).*Tlower;    % VECTOR OF LOWER TEMPERATURE LIMIT

% ***************************************************************
% SIMULATION
% ***************************************************************
for k = 1:n-1                % BEGINNING OF TIME LOOP
   
   if T(k) > Tupper          % CHECK FOR TURNING HEATER OFF
      Qf(k) = 0.0;           % SETS HEAT INPUT TO ZERO
   end                       % ENDS CHECK (OFF)
   if T(k) < Tlower          % CHECK FOR TURNING HEATER ON
      Qf(k) = Qfon;          % SETS HEAT INPUT TO Qfon
   end                       % ENDS CHECK (ON)
   if (T(k)<Tupper) & (T(k)>Tlower) & (k>2) % CHECK FOR INTERMEDIATE TEMP
      Qf(k) = Qf(k-1);       % KEEPS HEAT INPUT AT PREVIOUS HEAT INPUT
   end                       % ENDS CHECK (INTERMEDIATE TEMP)

% EULER INTEGRATION FOR DETERMINATION OF TEMPERATURE AND HEAT INPUT CHANGE
   
   if tauQ == 0.0            % CHECKS FOR FURNACE/ROOM HEAT DYNAMICS
      Qh(k) = Qf(k);         % IF NO DYNAMICS THEN ROOM HEAT EQUALS FURNACE
   else
      Qdot = (Qf(k) - Qh(k))/tauQ;
      Qh(k+1) = Qdot*delt + Qh(k);
   end                       % END IF (FURNACE/ROOM DYNAMICS)
   
   Tdot = (-(rho*Fair*Cv + UA)*(T(k) - Ta) + Qh(k))/(rho*V*Cv + Mwood*Cpwood);
   T(k+1) = Tdot*delt + T(k);
   
end                          % ENDS FOR (TIME LOOP)   

if tauQ == 0.0               % CHECKS IF UPDATING IS NECESSARY
   Qh(k+1) = Qh(k);          % PREDICTS HEAT INPUT FOR LAST ITERATION
end                          % END IF (UPDATE CHECK)

% ***************************************************************
% PLOTTING DATA
% ***************************************************************

subplot(2,1,1);plot(t,T,t,Tmax,'--',t,Tmin,'--');
axis([0 tend 15 25]);
title('Example 3.4');
ylabel('Room Temperature (degC)');xlabel('Time (Hours)');

subplot(2,1,2);plot(t,Qh);
axis([0 tend 0 2E6]);
ylabel('Heat Input To Room (cal/hr)');xlabel('Time (Hours)');

figure (1)
pause
close all