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/3/Le)*[7 -8 1;-8 16 -8;1 -8 7];   %elemental condutance matrix due to conduction
K_h=(ha*P*Le/30)*[4 2 -1;2 16 2;-1 2 4];   %elemental condutance matrix due to convection
K_e = K_C+K_h; % conductance matrix

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

%Contribution from concution side
for i=1:n
    K([2*i-1:2*i+1],[2*i-1:2*i+1]) =  K([2*i-1:2*i+1],[2*i-1:2*i+1])+K_e;
end

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

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

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

F = Fg+Fh;    %The RHS

%-------------solution------------------------%
K_act = K([2:2*n+1],[2:2*n+1]);                       %as the first node temperature is given
F_act = F([2:2*n+1])-K([2: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:2*n])];                      %Full temperature 
q1 = (K(1,[1:2*n+1])*T-F(1))/A;         

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

%plotting temperatures

s=0:0.05:1;                            %dimensionless co-ordinate
for i = 1:n                         %evaluation for element number
    for j=1:21                        %evaluate at each dimensionless co-ordinate
        x(20*(i-1)+j)=(i-1)*Le+Le*s(j);
        Temp(20*(i-1)+j)=T(2*i-1)*2*(s(j)-0.5)*(s(j)-1)+T(2*i)*4*s(j)*(1-s(j))+T(2*i+1)*2*s(j)*(s(j)-0.5);
    end
end

plot(x,Temp,'g');
xlabel('x (m)');
ylabel('Temperature (Deg C)');

hold on;

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

plot(x_m,T, 'o');

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:');
    
