clc;
clear all;

A = pi()/4*(0.5/12)^2;
P = pi()*(0.5/12);
L = 4/12;

k= 120;  
ha = 50;
hw = 100;
Ta = 72;
n = input('Number of elements:   ');

%--------------Calculations------------------------%

Le = L/n;                   %length of an element

K_C=(k*A/Le)*[1 -1;- 1 1];   %elemental condutance matrix due to conduction
K_h=(ha*P*Le/6)*[2 1;1 2];   %elemental condutance matrix due to convection
K_e = K_C+K_h; % conductance matrix

%-----------Assembly of Conductance matrix---------%
K = zeros(n+1,n+1);          %Global conductance matrix of DOF x DOF size

%Contribution from each element
for i=1:n
    K([i:i+1],[i:i+1]) =  K([i:i+1],[i:i+1])+K_e;
end

%updating the nth row, nth column by the convection term
K(n+1, n+1) =K(n+1, n+1)+A*hw;

%-----------gradient matrix------------------------%
Fg=zeros(n+1,1);
Fg(n+1) = 5.454;        %due to A*h_w*T_water

Fh_e = (ha*P*Ta*Le/2)*[1 1]';
Fh=zeros(n+1,1);
for i = 1:n
     Fh([i:i+1]) =  Fh([i:i+1])+Fh_e;
end 

F = Fg+Fh;    %The RHS

%-------------solution------------------------%
K_act = K([2:n+1],[2:n+1]);                       %as the first node temperature is given
F_act = F([2:n+1])-K([2:n+1],1)*180;     %active component of the gradient vector - note that given boundary condition on T is not zero
T_act=inv(K_act)*F_act;                          %solution for unknown temperatures

T = [180;T_act([1:n])];                      %Full temperature 
q1 = (K(1,[1:n+1])*T-F(1))/A;         

disp('    Node     Temperatures at nodes in Deg C');
noderesults = [single(1:n+1)',T];
disp(noderesults);
disp('Heat Flux in (W/m^2)');
disp(q1);

%plotting temperatures

x_m(1) = 0;
for i = 2:n+1
    x_m(i) = x_m(i-1)+Le;
end

plot(x_m,T, '-o');

xlabel('x (m)');
ylabel('Temperature (Deg C)');

hold on;

x_e(1) = 0;
for i = 1:101
    x_e(i) = (i-1)*L/100;
    T_e(i) = 72+0.76*exp(6.32*x_e(i))+107.24*exp(-6.32*x_e(i));
end

plot(x_e,T_e,'r:');
    
