clear; close all
%
%  Use 1-D thermal rod element for heat loss thru house wall
%
%  1/2" sheetrock, Kxx = 0.17 W/m*K
%  3.5" insulation, Kxx = 0.035 W/m*K
%  (neglect vapor barrier)
%  1" wood siding, Kxx = 0.12 W/m*K
%  Outside temp is 0C
%  Inside wall temp is 15C
%
nodes = [0  0.5  4   5]';  %% only x coord needed
nodes = nodes/39.37;        %% convert to m (sigh)
elems = [1 2
         2 3
         3 4];
nelems = size(elems,1);
Kxx = [0.17 0.035 0.12];
A = 1;      %% per m^2 of wall area
nnodes = length(nodes);
ndof = nnodes;
bigk = zeros(ndof,ndof);
for i=1:nelems
    n1 = elems(i,1);
    n2 = elems(i,2);
    kc = rodtk(nodes,n1,n2,Kxx(i),A);
    dof = [n1 n2];
    bigk(dof,dof) = bigk(dof,dof) + kc;
end
bc = [1  15         % node, temp
      4  0];
bcdof = bc(:,1);
freedof = 1:ndof;
freedof(bcdof) = [];
bigq = zeros(ndof,1);   %% no applied flux
bigt = zeros(ndof,1);
bigt(bcdof) = bc(:,2);
%     q_f             k_fb                    t_f
RHS = bigq(freedof) - bigk(freedof,bcdof)*bigt(bcdof);
bigt(freedof) = bigk(freedof,freedof)\RHS;
subplot(2,1,1)
plot(39.37*nodes,bigt,'-*');
ylabel('Temp, deg C')
subplot(2,1,2);
q = zeros(nelems,1);
for i=1:nelems
    n1 = elems(i,1);
    n2 = elems(i,2);
    L = nodes(n2) - nodes(n1);
    dTdx = (bigt(n2)-bigt(n1))/L;
    q(i) = -Kxx(i)*dTdx;        %% eqn 13.1.3
    plot(39.37*nodes([n1 n2],1),q(i)*[1 1],'-*');
    hold on
end
ylabel('Flux, watts/m^2');
xlabel('Position, in');
ylim(q(end)*[0.8 1.2]); %% expecting uniform flux
%
%   Cost calc
%
Awall = 20*8.5;                     %% Area of wall, ft^2
Awall = Awall * (12/39.37)^2;       %% m^2
watts = q(1)*Awall;                 %%
hrs = 24;                           %% calc for 24 hrs
kW = watts/1000;                    %% kW
kWhr = kW*hrs;                      %% kW-hrs
kWhr_cost = 0.10;                   %% 10 cents per hW-hr
day_cost = kWhr*kWhr_cost;
fprintf('Heating cost for this wall: $%.2f per day\n',day_cost);