% MorningStar Mechanical Room Model

%% Model Observations

%%   Initial Control Outputs

P1 = 0 ; % Solar Thermal

% Valves (0-10V Modulating control)
% Actual installation at MorningStar may not match schematic
% 0 corresponds to all flow from port B -> C
% 1 corresponds to all from from port A -> C
% 0-1 modulates mixed A & B flow percentages

H1 = 0; %Domestic hot water storage tank bypass, 0=total bypass, 1=all flow through tank
D1 = 1; %Heating loop storage tank bypass, 0=total bypass, 1=all flow through tank
R1 = .5; %Radiant zone 1 bypass, 0=total bypass, 1=all flow through loop
 %%  Properties of water at 40C (using system average)
rho = 1000; %kg/m^3
C_p = 4.179; %kJ/kgK
C = 60000; %This is a conversion factor to get from J/min to kW [L-s / m^3-min]
%%   MODEL PARAMETERS 
%Set Points
 T_hws = 60; %Deg C, Domestic Hot Water Supply, T_hws
 T_a_sp = 24.6; %Operating Temperature Set Point (Floor Temp can be 19-29C)
 T_rwsd = 1; %Half of anticipated radiant water temp drop to simulate average temperature
 
%Boundary Conditions
HighWater = 1; 

 
%Bypasses
 H1IO = 1; %Domestic water storage tank
 P1IO = 1; %Solar Thermal On/Off
 
%Excel Outputs
 xlsoutput = 0; %Final Day Array
 xlsoutput_e = 1; %Parametric Summary
 printscreen = 0; %Print all outputs to screen for final day
 printfigures = 1; %output figures for last day
 
 %Parameters
Cloud = 0;
Month = 3; %Jan=1 Mar=2 Apr=3 
F7 = 12; %L/min, Heating loop flow rate
 
%% Initial Conditions
 T6 = 30;    %C, starting entering water temp into solar collector
 T_Tank_u = 30; %C starting well-mixed tank temperature
 T_Tank_l = 30; %C starting well-mixed tank temperature

T_sb = 24;   %Slab Temps
T_s2 = 24;
T_s3 = 24;
T_s4 = 24;
T_s5 = 24;
T_st = 24;
T_a = 23;    %Air Temp
%% Loop Size and Time Step
tstep = 60;  %seconds, time step of each iteration, currently 1 min
n= 1440;     %loop size, 288 = 1 day @ tstep=300, 1440 = 1-day @ tstep=60
m=2;         %Day repeats to reach steady periodic, check convergence
b=20;        %Parametric Analysis Runs
%% Vector Pre-allocation

Plot_Q_t_bal = zeros(1,n);
Q_t_bal_sum = 0;
%Valves
e_i_d1 = zeros(1,n);
e_i_h1 = zeros(1,n);
e_i_r1 = zeros(1,n);

Plot_H1 = zeros(1,n);
Plot_D1 = zeros(1,n);
Plot_R1 = zeros(1,n);

%Temperatures

Plot_T1 = zeros(1,n);
Plot_T2 = zeros(1,n);
Plot_T3 = zeros(1,n);
Plot_T4 = zeros(1,n);
Plot_T5 = zeros(1,n);
Plot_T6 = zeros(1,n);
Plot_T7 = zeros(1,n);
Plot_T8 = zeros(1,n);
Plot_T9 = zeros(1,n);
Plot_T10 = zeros(1,n);
Plot_T11 = zeros(1,n);
Plot_T12 = zeros(1,n);
Plot_T13 = zeros(1,n);
Plot_T14 = zeros(1,n);
Plot_T15 = zeros(1,n);
Plot_T16 = zeros(1,n);

Plot_T_Tank = zeros(1,n);
Plot_T_Tank_l = zeros(1,n);
Plot_T_Tank_u = zeros(1,n);

%Flows

Plot_P1 = zeros(1,n);

Plot_F3 = zeros(1,n);
Plot_F6 = zeros(1,n);
Plot_F10 = zeros(1,n);

%Energy Calcs
Plot_E_tot = zeros(1,n);
Plot_E_input = zeros(1,n);
Plot_nu_st = zeros(1,n);
Plot_E_et = zeros(1,n);

Plot_E_mb = zeros(1,n);
Plot_E_wh = zeros(1,n);
Plot_E_shli = zeros(1,n);
Plot_E_whli = zeros(1,n);
Plot_q_1 = zeros(1,n);

Plot_i = zeros(1,n);
Plot_COP = zeros(1,n);
Plot_E_pv_total = zeros(1,n);
Plot_E_bal = zeros(1,n);

%Slab Model
Plot_T_st = zeros(1,n);
Plot_T_op = zeros(1,n);
Plot_dT_water = zeros(1,n);
Plot_T_a = zeros(1,n);
Plot_Q_ao = zeros(1,n);
Plot_T_slab = zeros(1,n);

%Excel outputs
T_op_avg_3  = zeros(1,b);
T_st_avg_3  = zeros(1,b);
T_a_avg_3  = zeros(1,b);
DailyBalance_3  = zeros(1,b);
E_mbt_3  = zeros(1,b);
E_wht_3  = zeros(1,b);
E_shl_3  = zeros(1,b);
E_whl_3  = zeros(1,b);
E_pv_total_3  = zeros(1,b);
E_et_total_3  = zeros(1,b);
COP_avg_3  = zeros(1,b);
TestVar = zeros(1,b);
T_Tank_avg_3 = zeros(1,b);
R1_avg_3 = zeros(1,b);

Plot_PMV = zeros(1,n);
Plot_PPD = zeros(1,n);

PPD_min_3= zeros(1,b);
PPD_max_3= zeros(1,b);
PMV_min_3= zeros(1,b);
PMV_max_3= zeros(1,b);
EPS_max_3= zeros(1,b);

%Storage Tank
Plot_Q_hx_l = zeros(1,n);
Plot_Q_hx_u = zeros(1,n);


%Zero summing variables
e_sum_d1=0;
e_sum_h1=0;
e_sum_r1=0;
E_tot=0;
E_whl=0;
E_shl=0;
m_w=0;
E_mbt=0;
E_wht=0;
E_pv1_total=0;
E_pv2_total=0;
E_et_tot=0;
q_1 = 0;

E_shl_2= 0;
E_whl_2 = 0;
E_demand_2 = 0;
E_mbt_2= 0;
E_wht_2 = 0;
E_tot_2= 0;
E_pv_total_2 = 0;
E_et_tot_2 = 0;
COP_avg_2 = 0;
DailyBalance_2 = 0;

T_op_sum = 0;
T_st_sum = 0;
T_a_sum = 0;
T_Tank_sum = 0;
R1_sum = 0;

bb = 0;
%% Read Heating Loads, Solar Irradiation, Ambient Temp

%February, coldest
if Month == 1

CL = 0.8; %Clothing Level for PMV calcs
T1 = 14.14; %C, based on (NREL, Hendron & Burch, 2007)
T_o = xlsread('MS_loads.xlsx','January','D2:D1441'); 

if HighWater == 1
F4 = xlsread('MS_loads.xlsx','January','C2:C1441'); 
else 
F4 = xlsread('MS_loads.xlsx','January','B2:B1441'); 
end

if Cloud == 1
I_s = xlsread('MS_loads.xlsx','January','G2:G1441'); %Includes tilt and clouds
else
I_s = xlsread('MS_loads.xlsx','January','F2:F1441'); %Sunny Day, Tilt, no cloud cover
end

else
end

%March
if Month == 2

CL = 0.7; %Clothing Level for PMV calcs
T1 = 14.43;
T_o = xlsread('MS_loads.xlsx','March','D2:D1441'); 

if HighWater == 1
F4 = xlsread('MS_loads.xlsx','March','C2:C1441'); 
else 
F4 = xlsread('MS_loads.xlsx','March','B2:B1441'); 
end

if Cloud == 1
I_s = xlsread('MS_loads.xlsx','March','G2:G1441'); %Includes tilt and clouds
else
I_s = xlsread('MS_loads.xlsx','March','F2:F1441'); %Sunny Day, Tilt, no cloud cover
end

else
end

%April, Mildest
if Month == 3

CL = 0.65; %Clothing Level for PMV calcs
T1 = 14.39;   
T_o = xlsread('MS_loads.xlsx','April','D2:D1441'); 

if HighWater == 1
F4 = xlsread('MS_loads.xlsx','April','C2:C1441'); 
else 
F4 = xlsread('MS_loads.xlsx','April','B2:B1441'); 
end

if Cloud == 1
I_s = xlsread('MS_loads.xlsx','April','G2:G1441'); %Includes tilt and clouds
else
I_s = xlsread('MS_loads.xlsx','April','F2:F1441'); %Sunny Day, Tilt, no cloud cover
end

else
end
%% Loops

F_stsp = 8; %L/min, Solar Thermal loop flow rate (Max set to 2GPM in specs)
T13 = 30; %Deg C, Radiant heating loop set point (54C Max)

%for F_stsp = [4,8,12]
%for T13 = [30,35,40,45,50]
%for F7 = [2-16]
for F7 = 16
    bb = 1 + bb;
for uu = 1:m % Number of days, repeating the same conditions. 
for i = 1:n %Going through n iterations
    
%% Domestic water sub-system PRE_STORAGE TANK
F1 = F4(i); %(L/min)Hot water demand is completely supplied from chilled water supply
F2 = F1*(1-H1);
F3 = F1*H1;
%% Radiant heating loop sub-system
%Start with micro-boiler set point and assumed flow rate

%F7 Assumed above
F6 = F7*D1  ;
F10 = F7*R1 ;
T15 = T13;

%If loops to avoid division by 0 when valves cut off flow

if F10 == 0
    T16 = T15;
else
    
    %% Floor Slab Model
%Model all three zones as a single zone. Not doing multiple iterations
%along the tube length due to the small temperature drop and multiple loops
%within the zone. 

Cp_slate = 760; %J/kg-K
rho_slate = 2700; %kg/m^3
k_slate = 2.01; %W/m-K
t_slate = 0.0254; %m (1")
n_slate = 5; %number of slices for analysis
R_slate = (t_slate/n_slate) / (k_slate); %m^2-C/W (resistance per linear meter along tubing)
U_slate = 1 / R_slate; %W/m^2-C
e_slate = .97;
sigma = 5.67*10^-8; %W/m2-K4
P_slate = 33.5; %m, Floor perimeter
Cp_air = 1005; %J/kg-K
rho_air = 1.205; %kg/m3
U_total = .3972; %W/m2-C
A_total = 175.2; %m2.
V_a = A_total * 3; %m3
D = (4 * A_total) / P_slate; %Hydraulic diameter
A = 58; %m^2, slab area
V_layer = A * (t_slate/n_slate); %m3, slice volume
Flow = F10; %L/min
rho_water = 1000; %kg/m^3
Cp_water = 4179; %J/kg-K
m_water = Flow * rho_water / 60000; %kg/s





%Assume air and walls are at same temperature
%Assume slices are at uniform temperature, given their small thicknesses


T_wi = T15-T_rwsd; %Drop X degrees to simulate an average temp because tube is modeled as single slice along tube direction
T_oa = T_o(i);

nn = 60;
tstepn = 1; %seconds

for ii = 1:nn
    
    h_r = e_slate * sigma * ((T_a+273)^2 + (T_st+273)^2) * ((T_a+273) + (T_st+273));
    h_c = (2.175/(D^.076))*((T_st-T_a)^.308);

    
    h_ta = h_r + h_c;
    hte = .5;
    
    
    Q_w1 = (m_water*Cp_water*(T_wi-T_sb))*hte;  
    Q_12 = U_slate * A * (T_sb - T_s2);
    Q_23 = U_slate * A * (T_s2 - T_s3);
    Q_34 = U_slate * A * (T_s3 - T_s4);
    Q_45 = U_slate * A * (T_s4 - T_s5);
    Q_5t = U_slate * A * (T_s5 - T_st);
    Q_ta = h_ta * A * (T_st - T_a);
    Q_ao = U_total * A_total * (T_a - T_oa);
    
    T_sb = T_sb + ((Q_w1-Q_12)*tstepn /(rho_slate*Cp_slate*V_layer)); 
    T_s2 = T_s2 + ((Q_12-Q_23)*tstepn /(rho_slate*Cp_slate*V_layer)); 
    T_s3 = T_s3 + ((Q_23-Q_34)*tstepn /(rho_slate*Cp_slate*V_layer)); 
    T_s4 = T_s4 + ((Q_34-Q_45)*tstepn /(rho_slate*Cp_slate*V_layer)); 
    T_s5 = T_s5 + ((Q_45-Q_5t)*tstepn /(rho_slate*Cp_slate*V_layer)); 
    T_st = T_st + ((Q_5t-Q_ta)*tstepn /(rho_slate*Cp_slate*V_layer)); 
    T_a = T_a + ((Q_ta-Q_ao) / (V_a * rho_air * Cp_air));
      
    T_wo = (T_wi - (Q_w1/(m_water*Cp_water)))+T_rwsd; %Outlet water temp
    dT_water = T_wi - (T_wo-T_rwsd); %Temperature drop across loop
    T_op = (T_st + T_a)/2; %Operating Temperature
end

    T16 = T_wo; 
    T_slab = (T_sb + T_s2 + T_s3 + T_s4 + T_s5 + T_st)/6;

end

    T9 =((T16*F10)+(T15*(F7-F10))) / F7  ;
%% Solar Thermal Array Model

% From the manufacturer:
% nu = 0.4916 - 1.9242 (W/m^2-C) * (Ti-Ta)/I
%Assume T6, then add energy based on assumed E_et (Solar input)

if P1 == 0
    F5 = 0;
else
    F5 = F_stsp; %L/min, Set Point
end


if I_s(i) == 0
    E_et=0;
    nu_st=0;
else 
nu_st = 0.4916-(1.9242*((T6-T_o(i))/I_s(i))) ; %Collector absorbtion efficiency

if nu_st < 0
    nu_st =0;
else 
end
E_et =(nu_st * 1.82 * I_s(i))/1000 ; % Absorption coefficient * Absorber Area [m^2] * Solar Irradiance (W/m^2)
end



if P1 == 1 %Is the pump running?
    T7 = T6 + ((C*E_et)/(C_p*F5*rho)); %Exiting water temp given solar energy input (Shoot for 200F or 93C)
else 
    T7=T6;
end
%% Storage Tank Heat Exchange

%Well-mixed sub-tanks with convective flows between

%Coil tube is composed of a steel tube with an enamel coating
%Radial heat exchange from the tube over two layers


o=60;
tstept = 1;

for ll = 1:o
    
D_i = .028; %m
D_o = .033; %m
t_steel = .002; %m (2 mm steel tubing)
t_coat = .0005; %m (0.5 mm Enamel Coating)
t_film = .001; %m (1 mm water film to simulate complex convective resistance)
r_water = D_i/2; %m

C_p_w = 4179; %kJ/kgK, water

k_steel = 45; %W/mk
k_coat = 0.8; %W/mk
k_water = 0.58; %W/m-K
j_u = F6 / (60 * pi * D_i / 4); %kg/s-m2
j_l = F5 / (60 * pi * D_i / 4); %kg/s-m2
mu = .0006532; %kg/m-s @ 40C

%Tank losses, assume inner surface temp = bulk tank section temp
k_ins = 0.020604; %W/mC
t_ins = 0.0762; %m
r_tank_o = .350; %m
r_tank_i = r_tank_o - t_ins;
h_t_u = .566; %m, tank section height
h_t_l = .97; %m, tank section height
R_ins_u = (log(r_tank_o/r_tank_i))/(2*pi*k_ins*h_t_u); %C/W
R_ins_l = (log(r_tank_o/r_tank_i))/(2*pi*k_ins*h_t_l); %C/W

if F6 == 0
    h_hxfluid_u = 1000000;
else
    h_hxfluid_u = 0.023 * ((k_water^(1-.33)*j_u^.8*C_p_w^.33)/(mu^(.8-.33)*D_i^.2));
end

if F5 == 0
    h_hxfluid_l = 1000000;
else
h_hxfluid_l = 0.023 * ((k_water^(1-.33)*j_l^.8*C_p_w^.33)/(mu^(.8-.33)*D_i^.2));
end

l_upper = 10.5; %m length of tubing, heating loop
l_lower = 14.45; %m length of tubing, solar 

A_s_upper = 1.1; %m^2, surface area
A_s_lower = 1.5; %m^2, surface area
A_s_in_u = pi*D_i*l_upper;
A_s_in_l = pi*D_i*l_lower;

%Total tank water mass = 305 kg, split two control volumes based on tank
%geometry
m_tank_storage = 305; %Total tank mass
m_tank_u = .37*m_tank_storage; %36.8% of mass is over upper coil
m_tank_l = .63*m_tank_storage;

 V_upper = 7.3; %liters
 V_lower = 10.1; %liters

T_hxdrop_l = 0; %Anticipated drop across coil to simulate bulk temp
T_hxdrop_u = 0;

F_conv = 0; %Combined convection flow
 %% Lower Tank Mass

%Lower Coil

R_hxfluid_l = 1 / (h_hxfluid_l * A_s_in_l); %Resistance between fluid and inner surface
R_steel_l = log((r_water+t_steel)/r_water)/(2*pi*k_steel*l_lower); %C/W
R_coat_l = log((r_water+t_steel+t_coat)/(r_water+t_steel))/(2*pi*k_coat*l_lower); %C/W
R_water_l = log((r_water+t_steel+t_coat+t_film)/(r_water+t_steel+t_coat))/(2*pi*k_water*l_lower); %C/W

R_tot_l = R_hxfluid_l + R_steel_l + R_coat_l + R_water_l;
U_tot_l = 1 / R_tot_l; %W/C

%Coil Heat Exchange
Q_hx_l = U_tot_l*((T7-T_hxdrop_l)-T_Tank_l); %Watts

%Exiting Water Temp
if F5 == 0
    Q_hx_l = 0;
    T6 = T7;
else
    T6 = T7 - (Q_hx_l / (F5/60 * C_p_w));
end

%Upper Water Draw, Incoming cold water from mains
Q_flow_l = (T_Tank_l - T9) * F6 * C_p_w/60; %Watts, due to water draw

%Standby Losses
Q_loss_l = (T_Tank_l - T_o(i))/R_ins_l; %Watts

%Convective exchange between tank sections
Q_conv_l = 0; %Watts

%Tank temperature balance, well mixed within section
T_Tank_l = T_Tank_l + (((Q_hx_l-Q_flow_l-Q_loss_l-Q_conv_l)*tstept) / (m_tank_l * C_p_w));

T8 = T_Tank_l;


%Average Tank Temp
T_Tank = T_Tank_l;

T2=T1;

T6_t = T6;
T7_t = T7;
T8_t = T8;
T9_t = T9;

%Storage Tank Energy Balance

if uu == m
 Q_t_in = ((T7-T6)*F5/60*C_p_w); %W
 Q_t_loss = (Q_loss_l);%W
 Q_t_out = ((T8-T9)*F6/60*C_p_w)+(Q_t_loss); %kW
Q_t_bal = Q_t_in - Q_t_out; %W
%Plot_Q_t_bal(i) = Q_t_bal;
Q_t_bal_sum = (Q_t_bal*tstept/3600000) + Q_t_bal_sum; %kWh
else
end

end
%% Micro-Boiler Model

T10 = ((T8*F6)+(T9*(F7-F6))) / F7     ;
T12 = T10   ;

if T10 > T13
    T14 = T10;
    E_mb = 0;
else
E_mb = ((T13-T10)*F7*rho*C_p) / C ; %Energy input into system, cost kW.
T14=T13;
end
%% Heat Pump Water Heater Model (Constant supply temp)

%T1 Assumed above, T2 solved in PRE-STORAGE section
if F1 == 0
    T3 = T1;
else
T3 = ((T1*F2)+(T2*F3)) / (F2+F3) ;
end


m_hpwh = 303; %kg or liters
%Mass of water incoming that needs to be heated:
m_hpwh_in = (F2+F3)*rho/C; %kg/s
COP_hpwh = 0.113*T_o(i) + 1.485; % From manufacturer

if T3>T_hws
    E_wh = 0;
    T4=T3;
else
E_wh = ((T_hws - T3)*m_hpwh_in*C_p)/COP_hpwh;
T4=T_hws;
end

T5=T4;
%% Solar PV Generation
% Array 1 = 12 panels
% Array 2 = 30 panels
% Panel Area = 1.25847 m^2 (1593x790x50mm)
% Total Area = 52.855 m^2
% BP Solar 7190 panel efficiency = 15.1%

%Array 1 is connected to inverter 1, which has a capacity of 1800W
%Array 2 is connected to inverter 2, which has a capacity of 5700W

E_pv1 = (.151 * 15.10164 * I_s(i))/1000; %kW generation, array 1

if E_pv1 > 1.800 %Do not exceed inverter capacity
    E_pv1 = 1.800;
else
end 

E_pv2 = (.151 * 37.7541 * I_s(i))/1000; %kW generation, array 2

if E_pv2 > 5.700 %Do not exceed inverter capacity
    E_pv2 = 5.700;
else
end 

E_pv = E_pv1 + E_pv2;
E_pv1_total = (E_pv1 * (tstep/3600)) + E_pv1_total; %kWh Total PV Generation
E_pv2_total = (E_pv2 * (tstep/3600)) + E_pv2_total; %kWh Total PV Generation
E_pv_total = E_pv1_total + E_pv2_total;
E_pvi_total = E_pv * (tstep/3600);
%% Total Energy Usage

E_input = E_mb + E_wh; %Micro boiler and HPWH load
E_tot = (E_input * (tstep/3600)) + E_tot ; %kW*iterations
m_w = (F4(i)*tstep/60) + m_w; %Total water use

E_mbt = (E_mb * (tstep/3600)) + E_mbt; %Total micro-boiler energy input
E_wht = (E_wh * (tstep/3600)) + E_wht; %Total HPWH energy input

E_shl = ((Q_ao/1000)*(tstep/3600))+ E_shl ; %Space heating load
E_whl = ((F4(i)*tstep/60)*(rho/1000)*(T_hws-T1)* C_p/3600)+ E_whl ; %Water heating load %kJ*3600=kWh
E_demand = E_shl + E_whl;

COP_avg = E_demand / E_tot;

%Instantaneous Loads

E_mbi = (E_mb * (tstep/3600)) ; %Total micro-boiler energy input
E_whi = (E_wh * (tstep/3600)) ; %Total HPWH energy input

E_shli = ((Q_ao/1000)*(tstep/3600)) ; %Space heating load
E_whli = ((F4(i)*tstep/60)*(rho/1000)*(T_hws-T1)* C_p/3600) ; %Water heating load %kJ*3600=kWh

COP = (E_shli + E_whli) / (E_mbi + E_whi);
DailyBalance = E_tot - E_pv_total; %kWh

%Energy balance with grid

E_bal = E_mbi + E_whi - E_pvi_total; %kWh

%Solar Thermal Generation

E_et_tot = (E_et * (tstep/3600))+E_et_tot; %kWh

%% PMV and PPD Calculator


%Inputs

CLO = CL;    %Clothing level
MET = 1.1;    %Metabolic rate
WME = 0;    %External Work
TA = T_a;   %Air Temperature
TR = T_st;  %Mean Radiant Temperature, Assume top of slab temp
VEL = 0.15;   %Relative Air Velocity
RH = 50;    %Relative Humidity

% Definitions

FNPS = exp(16.6536-4030.183/(TA+235));  %Saturated Vapor Pressure, KPa
PA = RH*10*FNPS;                        %Vapor Pressure, Pa
ICL = .155 * CLO;                       %Thermal Insulation of Clothing, m2K/W
M = MET * 58.15;                        %Metabolic Rate, W/m2
W = WME * 58.15;                        %External Work, W/m2
MW = M+W;                               %Internal Heat Production

if ICL < .078                           %Clothing Area Factor
    FCL = 1+1.29*ICL;
else
    FCL = 1.05+.645*ICL;
end

HCF = 12.1*VEL^.5;                       %h_conv, forced convection
TAA = TA+273;                            %Air Temp, K
TRA = TR+273;                            %Radiant Temp, K


% Surface Temp of clothing, by iteration

TCLA = TAA + (35.5-TA) / (3.5*(6.45*ICL+.1));   %First guess

P1 = ICL * FCL;                                 %Calculation Terms
P2 = P1 * 3.96;
P3 = P1 * 100;
P4 = P1* TAA;
P5 = 308.7 - .028*MW+P2*(TRA/100)^4;
XN = TCLA/100;
XF = TCLA/50;                                   %Pass initial stop criteria


N=0;                                            %Number of iterations
EPS = .00015;                                   %Stop criteria

while abs(XN-XF) > EPS            
XF = (XF+XN)/2;

HCN = 2.38*abs(100*XF-TAA)^.25;                 %h_c, natural convection

if HCF > HCN
    HC = HCF;
else
    HC = HCN;
end

XN = (P5+P4*HC-P2*XF^4)/(100+P3*HC);
N=N+1;
end

TCL = 100*XN-273;                                %Surface Temp of clothing

% Heat Loss Through Clothing

HL1 = 3.05*.001*(5733-6.99*MW-PA);      %Heat loss diff. through skin
if MW >58.15
    HL2 = .42*(MW-58.15);               %Heat loss by sweating
else
    HL2=0;
end
HL3 = 1.7 * .00001 * M * (5867-PA);     %Latent respiration heat loss
HL4 = .0014 * M * (34-TA);              %Dry respiration heat loss
HL5 = 3.96*FCL*(XN^4-(TRA/100)^4);      %Heat loss by radiation
HL6 = FCL * HC * (TCL-TA);              %Convective heat loss

% Calculate PMV and PPD

TS = .303 * exp(-.036*M)+.028;              %Thermal sensation trans. coeff.
PMV = TS * (MW-HL1-HL2-HL3-HL4-HL5-HL6);    %Predicted Mean Vote
PPD = 100-95*exp(-.03353*PMV^4-.2179*PMV^2);%Predicted Percent Dissatisfied
%% Final Simulation Day Summations

if uu ==m

    
T_op_sum = T_op + T_op_sum;
T_op_avg = T_op_sum / (i);
T_st_sum = T_st + T_st_sum;
T_st_avg = T_st_sum / (i);
T_a_sum = T_a + T_a_sum;
T_a_avg = T_a_sum / (i);
T_Tank_sum = T_Tank + T_Tank_sum;
T_Tank_avg = T_Tank_sum / (i);
R1_sum = R1_sum + R1;
R1_avg = R1_sum/(i);

E_shl_2 = ((Q_ao/1000)*(tstep/3600))+ E_shl_2; %Space Heating Load
E_whl_2 = ((F4(i)*tstep/60)*(rho/1000)*(T_hws-T1)* C_p/3600)+ E_whl_2; %Water Heating Load 
E_demand_2 = E_shl_2+E_whl_2; %Total Energy Demand
E_mbt_2 = (E_mb * (tstep/3600)) + E_mbt_2; %Micro Boiler Usage
E_wht_2 = (E_wh * (tstep/3600)) + E_wht_2; %Water Heater Usage
E_tot_2 = (E_input * (tstep/3600)) + E_tot_2; %Total Energy Usage
E_pv_total_2 = ((E_pv1 + E_pv2)*(tstep/3600)) + E_pv_total_2; %Total PV Generation
E_et_tot_2 = (E_et * (tstep/3600))+E_et_tot_2; %Solar Thermal Generation
COP_avg_2 = E_demand_2 / E_tot_2;
DailyBalance_2 = E_tot_2 - E_pv_total_2; %Net energy into system    
    

%% Outputs variables for parametric analysis


T_op_avg_3(bb) = T_op_avg;
T_st_avg_3(bb) = T_st_avg; 
T_a_avg_3(bb) = T_a_avg;
DailyBalance_3(bb) = DailyBalance_2;
E_mbt_3(bb)=E_mbt_2;
E_wht_3(bb)=E_wht_2;
E_shl_3(bb)=E_shl_2;
E_whl_3(bb)=E_whl_2;
E_pv_total_3(bb)=E_pv_total_2;
E_et_total_3(bb)=E_et_tot_2;
COP_avg_3(bb)=COP_avg_2;
TestVar(bb) = F7;
T_Tank_avg_3(bb) = T_Tank_avg;
R1_avg_3(bb) = R1_avg;

else
end
%% Output results


if uu==m
 if printscreen == 1   
%Iteration number
fprintf(1,'Iteration number ') ;
fprintf(1,'%3.2f \n', i) ;

%Floop loop temps in and out
fprintf(1,'  \n') ;
fprintf(1,'Temperatures beginning and ending floor loop: \n') ;
fprintf(1,'%3.2f', T15) ;
fprintf(1,'C T15, Radiant Floor In \n') ;
fprintf(1,'%3.2f', T16) ;
fprintf(1,'C T16, Radiant Floor Out \n') ;
fprintf(1,'%3.2f', T9) ;
fprintf(1,'C T9, Mixed end, storage tank in \n') ;

%Floor Loop
fprintf(1,'  \n') ;
fprintf(1,'Temperatures drops across radiant floor loops: \n') ;
fprintf(1,'%3.2f', dT_water) ;
fprintf(1,'C dT_water \n') ;


%Domestic Water Temps
fprintf(1,'  \n') ;
fprintf(1,'Domestic Water Temps: \n') ;
fprintf(1,'%3.2f', T1) ;
fprintf(1,'C T1p \n') ;
fprintf(1,'%3.2f', T2) ;
fprintf(1,'C T2 \n') ;
fprintf(1,'%3.2f', T3) ;
fprintf(1,'C T3 \n') ;
fprintf(1,'%3.2f', T4) ;
fprintf(1,'C T4 \n') ;

%Heating Loop Temps
fprintf(1,'  \n') ;
fprintf(1,'Heating Loop Temps: \n') ;
fprintf(1,'%3.2f', T9) ;
fprintf(1,'C T9 \n') ;
fprintf(1,'%3.2f', T8) ;
fprintf(1,'C T8 \n') ;
fprintf(1,'%3.2f', T10) ;
fprintf(1,'C T10 \n') ;
fprintf(1,'%3.2f', T13) ;
fprintf(1,'C T13 \n') ;

%Solar Loop Temps
fprintf(1,'  \n') ;
fprintf(1,'Solar Loop Temps: \n') ;
fprintf(1,'%3.2f', T7) ;
fprintf(1,'C T7 \n') ;
fprintf(1,'%3.2f', T6) ;
fprintf(1,'C T6 \n') ;

%Storage Tank Temps
fprintf(1,'  \n') ;
fprintf(1,'Storage Tank Temps: \n') ;
fprintf(1,'%3.2f', T1) ;
fprintf(1,'C T1 \n') ;
fprintf(1,'%3.2f', T_Tank_l) ;
fprintf(1,'C Lower Tank \n') ;
fprintf(1,'%3.2f', T_Tank_u) ;
fprintf(1,'C Upper Tank \n') ;
fprintf(1,'%3.2f', T2) ;
fprintf(1,'C T2 \n') ;
    
%Total Energy Use
fprintf(1,'  \n') ;
fprintf(1,'Space Heating Load: \n') ;
fprintf(1,'%3.5f', E_shl_2) ;
fprintf(1,' kWh  \n') ;
fprintf(1,'Water Heating Load: \n') ;
fprintf(1,'%3.5f', E_whl_2) ;
fprintf(1,' kWh  \n') ;
fprintf(1,'Total Energy Demand: \n') ;
fprintf(1,'%3.5f', E_demand_2) ;
fprintf(1,' kWh  \n') ;
fprintf(1,'Micro-Boiler Enery Cost: \n') ;
fprintf(1,'%3.5f', E_mbt_2) ;
fprintf(1,' kWh  \n') ;
fprintf(1,'HPWH Enery Cost: \n') ;
fprintf(1,'%3.5f', E_wht_2) ;
fprintf(1,' kWh  \n') ;
fprintf(1,'Total Enery Cost: \n') ;
fprintf(1,'%3.5f', E_tot_2) ;
fprintf(1,' kWh  \n') ;
fprintf(1,'Total PV Generation: \n') ;
fprintf(1,'%3.5f', E_pv_total_2) ;
fprintf(1,' kWh  \n') ;
fprintf(1,'Total Solar Thermal Generation: \n') ;
fprintf(1,'%3.5f', E_et_tot_2) ;
fprintf(1,' kWh  \n') ;
fprintf(1,'Coefficient of Performance: \n') ;
fprintf(1,'%3.5f', COP_avg_2) ;
fprintf(1,'  \n') ;
fprintf(1,'Net Daily Energy Balance: \n') ;
fprintf(1,'%3.5f', DailyBalance_2) ;
fprintf(1,'kWh  \n') ;


%Control Outputs
fprintf(1,'  \n') ;
fprintf(1,'Control Outputs: \n') ;
fprintf(1,'%3.2f', H1) ;
fprintf(1,' H1, valve position going into HPWH \n') ;
fprintf(1,'%3.2f', D1) ;
fprintf(1,' D1, valve position for heating loop to bypass storage tank \n') ;
fprintf(1,'%3.2f', R1) ;
fprintf(1,' R1, valve position modulating flow radiant floor \n') ;
fprintf(1,'%3.2f', P1) ;
fprintf(1,' P1, Pump status: solar hot water collector \n') ;
fprintf(1,'  \n') ;
 else 
 end
 
else
end
%% Update control valve positions
% Update assumed initial temps (T14, T6)
% Update pump status

%Implement PID control
% u(t) = MV(t) = K_p * e(t) + K_i * Int_0_t(e(tau)dtau) + K_d * d/dt e(t)
% K_i = Integral Gain
% K_p = Proportional Gain
% K_d = Derivative Gain
% e = error = Set Point - Present Value
% t = instantaneous time = iteration number
% tau = Variable of integration 0 < tau < t
%% Valve H1 Pre-HPWH, after storage tank in domestic supply

H1=0;

%% Valve D1, Heating Loop Storage Tank Bypass

D1=1;

%% Valve R1, Radiant Floor Flow Modulation
% Trying to achieve the design deltaT of 10F or 5.6C

   
%  K_p_r1 = (.089285/2);      %Gain, Proportional term
 K_p_r1 = .09;      %Gain, Proportional term
 K_i_r1 = 0.003;      %Gain, Integral term
%  K_d_r1 = 0;      %Gain, Derivative term
 er_r1 = T_a_sp - T_op;   %Current error term
 
 
 for u = 1:n %Build memory of error terms for integral and derivative calculations
    e_i_r1(u) = er_r1 ;
 end
 
%  if T15 <= (T16+10)
%     R1 = 1;
%  else

    %PI Control
    e_sum_r1 = e_sum_r1 + e_i_r1(u); %Sum of error terms for integral term
    R1 = .5 + (K_p_r1 * er_r1) + (K_i_r1 * e_sum_r1); % - (K_d_r1 * (e_i_r1(u)-e_i_r1(u-1)));
    
%  end

    if R1>1
    R1=1;
    else
    end
    
     if R1<0.1
     R1=0;
    else
     end  
%% Pump P1 Solar Thermal
% Should be turned off when there is no longer and thermal input, ie. E_et
% = 0, OR when the collector outlet temp is getting too hot, 
% and cannot take any more heat.
if P1IO == 1

if E_et > 0 && T7 < 110
    P1 = 1;
else 
    P1 = 0;
    T6 = T_Tank_l;
    T7 = T_Tank_l;
end

else
    P1=0;
end
%% Collect Data for plotting

if uu==m
    
Plot_H1(i) = H1;
Plot_D1(i) = D1;
Plot_R1(i) = R1;
Plot_F3(i) = F3;
Plot_P1(i) = P1;

Plot_T1(i) = T1;
Plot_T2(i) = T2;
Plot_T3(i) = T3;
Plot_T4(i) = T4;
Plot_T6(i) = T6;
Plot_T7(i) = T7;
Plot_T8(i) = T8;
Plot_T9(i) = T9;
Plot_T10(i) = T10;
Plot_T13(i) = T13;
Plot_T14(i) = T14;
Plot_T15(i) = T15;
Plot_T16(i) = T16;

Plot_T_Tank(i) = T_Tank;
Plot_T_Tank_l(i) = T_Tank_l;
Plot_T_Tank_u(i) = T_Tank_u;

Plot_F6(i) = F6;
Plot_F10(i) = F10;

Plot_E_tot(i) = E_tot;
Plot_E_input(i) = E_input;
Plot_nu_st(i) = nu_st;
Plot_E_et(i) = E_et;

Plot_E_mb(i)=E_mb;
Plot_E_wh(i)=E_wh;
Plot_E_shli(i)=E_shli;
Plot_E_whli(i)=E_whli;
Plot_q_1(i)=q_1;

Plot_i(i)=i;
Plot_COP(i)=COP;
Plot_E_pv_total(i)=E_pv;
Plot_E_bal(i) = E_bal;

Plot_T_op(i) = T_op;
Plot_T_st(i) = T_st;
Plot_dT_water(i) = dT_water;
Plot_T_a(i) = T_a;
Plot_Q_ao(i) = Q_ao/1000;
Plot_T_slab(i)=T_slab;

Plot_PMV(i) = PMV;
Plot_PPD(i) = PPD;

PPD_min_3(bb) = min(Plot_PPD);
PPD_max_3(bb) = max(Plot_PPD); %Must be <10%
PMV_min_3(bb) = min(Plot_PMV); %No less than 0.5
PMV_max_3(bb) = max(Plot_PMV); %No more than 0.5


Plot_Q_hx_l(i) = Q_hx_l/1000; %W to kW
Plot_Q_hx_u(i) = 0;


%Tests for convergence
EPS_T_Tank_l=(Plot_T_Tank_l(n)-Plot_T_Tank_l(1))/(Plot_T_Tank_l(n)); %Storage Tank
EPS_T_Tank_u=(Plot_T_Tank_u(n)-Plot_T_Tank_u(1))/(Plot_T_Tank_u(n)); %Storage Tank
EPS_T_a=(Plot_T_a(n)-Plot_T_a(1))/(Plot_T_a(n)); %Room Air Temp
EPS_T_st=(Plot_T_st(n)-Plot_T_st(1))/(Plot_T_st(n)); %Slab Temp
EPS_E_mb=(Plot_E_mb(n)-Plot_E_mb(1))/(Plot_E_mb(n)); %Micro-boiler Energy Consumption
EPS = [abs(EPS_E_mb),abs(EPS_T_st),abs(EPS_T_a),abs(EPS_T_Tank_l),abs(EPS_T_Tank_u)];
EPS_max=max(EPS);

EPS_max_3(bb) = EPS_max;

else
end
%%
end % day loop
end % Day repeat loop
%% Experimental Results Print to Screen

fprintf(1,'  \n') ;
fprintf(1,'Iteration: \n') ;
fprintf(1,'%3.5f', bb );
fprintf(1,'  \n') ;
fprintf(1,'Space Heating Load: \n') ;
fprintf(1,'%3.5f', E_shl_3(bb)) ;
fprintf(1,' kWh  \n') ;
fprintf(1,'Water Heating Load: \n') ;
fprintf(1,'%3.5f', E_whl_3(bb)) ;
fprintf(1,' kWh  \n') ;
fprintf(1,'Micro-Boiler Enery Cost: \n') ;
fprintf(1,'%3.5f', E_mbt_3(bb)) ;
fprintf(1,' kWh  \n') ;
fprintf(1,'HPWH Enery Cost: \n') ;
fprintf(1,'%3.5f', E_wht_3(bb)) ;
fprintf(1,' kWh  \n') ;
fprintf(1,'Total PV Generation: \n') ;
fprintf(1,'%3.5f', E_pv_total_3(bb)) ;
fprintf(1,' kWh  \n') ;
fprintf(1,'Total Solar Thermal Generation: \n') ;
fprintf(1,'%3.5f', E_et_total_3(bb)) ;
fprintf(1,' kWh  \n') ;
fprintf(1,'Coefficient of Performance: \n') ;
fprintf(1,'%3.5f', COP_avg_3(bb)) ;
fprintf(1,'  \n') ;
fprintf(1,'Net Daily Energy Balance: \n') ;
fprintf(1,'%3.5f', DailyBalance_3(bb)) ;
fprintf(1,'kWh  \n') ;
fprintf(1,'PMV Min: \n') ;
fprintf(1,'%3.5f', PMV_min_3(bb)) ;
fprintf(1,'  \n') ;
fprintf(1,'PMV Max: \n') ;
fprintf(1,'%3.5f', PMV_max_3(bb)) ;
fprintf(1,'  \n') ;
fprintf(1,'PPD min: \n') ;
fprintf(1,'%3.5f', PPD_min_3(bb)) ;
fprintf(1,'  \n') ;
fprintf(1,'PPD max: \n') ;
fprintf(1,'%3.5f', PPD_max_3(bb)) ;
fprintf(1,'  \n') ;
fprintf(1,'R1 Average \n') ;
fprintf(1,'%3.5f', R1_avg_3(bb)) ;
fprintf(1,'  \n') ;
fprintf(1,'Storage Tank Balance \n') ;
fprintf(1,'%3.5f', Q_t_bal_sum) ;
fprintf(1,'  \n') ;
fprintf(1,'Convergence Checks: \n') ;
fprintf(1,'%3.5f', EPS_T_Tank_l) ;
fprintf(1,'  \n') ;
 fprintf(1,'%3.5f', EPS_T_Tank_u) ;
 fprintf(1,'  \n') ;
 fprintf(1,'%3.5f', EPS_T_a) ;
 fprintf(1,'  \n') ;
 fprintf(1,'%3.5f', EPS_T_st) ;
 fprintf(1,'  \n') ;
 fprintf(1,'%3.5f', EPS_E_mb) ;
 fprintf(1,'  \n') ;
%% Zero Summing Variables for next parametric iteration

T_op_sum = 0;
T_op_avg = 0;
T_st_sum = 0;
T_st_avg = 0;
T_a_sum = 0;
T_a_avg = 0;
T_Tank_sum = 0;
T_Tank_avg = 0;
R1_avg = 0;
R1_sum = 0;

E_shl_2 = 0; %Space Heating Load
E_whl_2 = 0; %Water Heating Load 
E_demand_2 = 0; %Total Energy Demand
E_mbt_2 = 0; %Micro Boiler Usage
E_wht_2 = 0; %Water Heater Usage
E_tot_2 = 0; %Total Energy Usage
E_pv_total_2 = 0; %Total PV Generation
E_et_tot_2 = 0; %Solar Thermal Generation
COP_avg_2 = 0;
DailyBalance_2 = 0; %Net energy into system
Q_t_bal_sum = 0;

%%
end % Experiment Loop
%% Experiment Results to Excel

if xlsoutput_e==1

xlswrite('MSoutput2.xlsx',{'Heating Loop Flow Rate'},'A1:A1');
xlswrite('MSoutput2.xlsx',TestVar','A2:A21'); 
%Metrics
xlswrite('MSoutput2.xlsx',{'Oper. Temp, avg'},'B1:B1');
xlswrite('MSoutput2.xlsx',T_op_avg_3','B2:B21');

xlswrite('MSoutput2.xlsx',{'Slab Temp, avg'},'C1:C1');
xlswrite('MSoutput2.xlsx',T_st_avg_3','C2:C21');

xlswrite('MSoutput2.xlsx',{'Room Temp, avg'},'D1:D1');
xlswrite('MSoutput2.xlsx',T_a_avg_3','D2:D21');

xlswrite('MSoutput2.xlsx',{'Daily Energy Bal'},'E1:E1');
xlswrite('MSoutput2.xlsx',DailyBalance_3','E2:E21');

xlswrite('MSoutput2.xlsx',{'Micro-boiler Energy'},'F1:F1');
xlswrite('MSoutput2.xlsx',E_mbt_3','F2:F21');

xlswrite('MSoutput2.xlsx',{'HPWH Energy'},'G1:G1');
xlswrite('MSoutput2.xlsx',E_wht_3','G2:G21');

xlswrite('MSoutput2.xlsx',{'Space Loads'},'H1:H1');
xlswrite('MSoutput2.xlsx',E_shl_3','H2:H21');

xlswrite('MSoutput2.xlsx',{'Water Heating Loads'},'I1:I1');
xlswrite('MSoutput2.xlsx',E_whl_3','I2:I21');

xlswrite('MSoutput2.xlsx',{'PV Generation'},'J1:J1');
xlswrite('MSoutput2.xlsx',E_pv_total_3','J2:J21');

xlswrite('MSoutput2.xlsx',{'Solar Thermal Generation'},'K1:K1');
xlswrite('MSoutput2.xlsx',E_et_total_3','K2:K21');

xlswrite('MSoutput2.xlsx',{'COP'},'L1:L1');
xlswrite('MSoutput2.xlsx',COP_avg_3','L2:L21');

xlswrite('MSoutput2.xlsx',{'PPD max'},'M1:M1');
xlswrite('MSoutput2.xlsx',PPD_max_3','M2:M21');

xlswrite('MSoutput2.xlsx',{'Tank Temp, Avg'},'N1:N1');
xlswrite('MSoutput2.xlsx',T_Tank_avg_3','N2:N21');

xlswrite('MSoutput2.xlsx',{'Valve R1, Avg'},'O1:O1');
xlswrite('MSoutput2.xlsx',R1_avg_3','O2:O21');

xlswrite('MSoutput2.xlsx',{'Convergence'},'P1:P1');
xlswrite('MSoutput2.xlsx',EPS_max_3','P2:P21');

xlswrite('MSoutput2.xlsx',{'Tank Balance'},'Q1:Q1');
xlswrite('MSoutput2.xlsx',Q_t_bal_sum','Q2:Q21');


%Boundary Conditions
xlswrite('MSoutput2.xlsx',{'Cloudy?'},'S2:S2');
xlswrite('MSoutput2.xlsx',Cloud,'T2:T2');

xlswrite('MSoutput2.xlsx',{'Simulated Day'},'S4:S4');
xlswrite('MSoutput2.xlsx',Month,'T4:T4');

%Modes
xlswrite('MSoutput2.xlsx',{'Domestic Water Storage Bypass'},'S5:S5');
xlswrite('MSoutput2.xlsx',H1IO,'T5:T5');


%Set Points
xlswrite('MSoutput2.xlsx',{'Micro-boiler Set Point'},'S7:S7');
xlswrite('MSoutput2.xlsx',T13,'T7:T7');

xlswrite('MSoutput2.xlsx',{'Operating Temp Set Point'},'S8:S8');
xlswrite('MSoutput2.xlsx',T_a_sp,'T8:T8');

xlswrite('MSoutput2.xlsx',{'Solar Thermal Flow Rate'},'S9:S9');
xlswrite('MSoutput2.xlsx',F_stsp,'T9:T9');

xlswrite('MSoutput2.xlsx',{'Heating Loop Flow Rate'},'S10:S10');
xlswrite('MSoutput2.xlsx',F7,'T10:T10');

else
end
%
%% Output Figures 
if printfigures == 1
    
figure(1)
 plot(Plot_PPD)
 title('Percent People Dissatisfied');
figure(2)
 plot(Plot_T_op)
  title('Operating Temp');
figure(3)
 plot(Plot_D1)
  title('Valve D1 Control Signal');
figure(4)
 plot(Plot_PMV)
  title('PMV');
figure(5)
 plot(Plot_nu_st)
  title('Solar Thermal Efficiency');
% figure(6)
%  plot(Plot_T_a)
%   title('Room Air Temperature');
  figure(6)
 plot(Plot_T9)
  title('T9');
figure(7)
 plot(Plot_T_Tank_u)
  title('Upper Storage Tank Temperature');
figure(8)
 plot(Plot_E_mb)
  title('Micro-Boiler Energy Use');
figure(9)
 plot(Plot_R1)
  title('Valve R1 Control Signal');
% figure(10)
%  plot(Plot_T_st)
%   title('Slab Surface Temperature');
figure(10)
 plot(Plot_T7)
  title('Lower Coil Inlet');
figure(11)
 plot(Plot_E_wh)
  title('HPWH Energy Use');
figure(12)
 plot(Plot_E_et)
  title('Solar Thermal Generation');
figure(13)
 plot(Plot_Q_hx_l)
  title('Lower Loop Heat Exchange');
figure(14)
 plot(Plot_T_Tank_l)
  title('Lower Tank Temp');
figure(16)
 plot(Plot_F3)
  title('Water Draw From Storage Tank');
else
end
%
%% Output final day loop to Excel
if xlsoutput == 1 
    
xlswrite('MSoutput.xlsx',{'Set Points'},'A1:A1');
xlswrite('MSoutput.xlsx',{'Micro Boiler Set Point'},'A2:A2');
xlswrite('MSoutput.xlsx',T13,'A3:A3');
xlswrite('MSoutput.xlsx',{'Radiant Loop Flow Rate'},'A4:A4');
xlswrite('MSoutput.xlsx',F7,'A5:A5');
xlswrite('MSoutput.xlsx',{'Domestic Hot Water Supply Temp'},'A6:A6');
xlswrite('MSoutput.xlsx',T_hws,'A7:A7');
xlswrite('MSoutput.xlsx',{'Domestic Chilled Water Supply Temp'},'A8:A8');
xlswrite('MSoutput.xlsx',T1,'A9:A9');
xlswrite('MSoutput.xlsx',{'Solar Thermal Flow Rate'},'A10:A10');
xlswrite('MSoutput.xlsx',F_stsp,'A11:A11');
xlswrite('MSoutput.xlsx',{'Operating Temperature Set Point'},'A12:A12');
xlswrite('MSoutput.xlsx',T_a_sp,'A13:A13');
xlswrite('MSoutput.xlsx',{'Radiant loop valve controlled by operating temperature'},'A14:A14');
  
xlswrite('MSoutput.xlsx',{'Iteration'},'B1:B1');
xlswrite('MSoutput.xlsx',Plot_i','B2:B1441');

xlswrite('MSoutput.xlsx',{'HPWH Energy Input'},'C1:C1');
xlswrite('MSoutput.xlsx',Plot_E_wh','C2:C1441');

xlswrite('MSoutput.xlsx',{'Space Heating Demand'},'D1:D1');
xlswrite('MSoutput.xlsx',Plot_E_shli','D2:D1441');

xlswrite('MSoutput.xlsx',{'Water Heating Demand'},'E1:E1');
xlswrite('MSoutput.xlsx',Plot_E_whli','E2:E1441');

xlswrite('MSoutput.xlsx',{'Solar Thermal Array Generation'},'F1:F1');
xlswrite('MSoutput.xlsx',Plot_E_et','F2:F1441');

xlswrite('MSoutput.xlsx',{'Storage Tank Temperature'},'G1:G1');
xlswrite('MSoutput.xlsx',Plot_T_Tank','G2:G1441');

xlswrite('MSoutput.xlsx',{'PV Array Generation'},'H1:H1');
xlswrite('MSoutput.xlsx',Plot_E_pv_total','H2:H1441');

xlswrite('MSoutput.xlsx',{'Energy Balance'},'I1:I1');
xlswrite('MSoutput.xlsx',Plot_E_bal','I2:I1441');

xlswrite('MSoutput.xlsx',{'Storage Tank Energy Balance'},'J1:J1');
xlswrite('MSoutput.xlsx',Plot_Q_t_bal','J2:J1441');

xlswrite('MSoutput.xlsx',{'Radiant Floor Loop Output'},'K1:K1');
xlswrite('MSoutput.xlsx',Plot_T16','K2:K1441');

xlswrite('MSoutput.xlsx',{'Valve D1'},'L1:L1');
xlswrite('MSoutput.xlsx',Plot_D1','L2:L1441');

xlswrite('MSoutput.xlsx',{'Instantaneous COP'},'M1:M1');
xlswrite('MSoutput.xlsx',Plot_COP','M2:M1441');

xlswrite('MSoutput.xlsx',{'Micro-boiler Energy Input'},'N1:N1');
xlswrite('MSoutput.xlsx',Plot_E_mb','N2:N1441');

xlswrite('MSoutput.xlsx',{'Radiant Floor Flow Rate'},'O1:O1');
xlswrite('MSoutput.xlsx',Plot_F10','O2:O1441');

xlswrite('MSoutput.xlsx',{'Operating Temperature'},'P1:P1');
xlswrite('MSoutput.xlsx',Plot_T_op','P2:P1441');

xlswrite('MSoutput.xlsx',{'Valve R1'},'Q1:Q1');
xlswrite('MSoutput.xlsx',Plot_R1','Q2:Q1441');

xlswrite('MSoutput.xlsx',{'External Load'},'R1:R1');
xlswrite('MSoutput.xlsx',Plot_Q_ao','R2:R1441');
else
end 
%
%% 
fprintf(1,'done!');
