clear all;   clc;

T=[]; % Temperature storage for all unknowns i.e. n-1 nodes
T_all=[]; % Including the known temperature
a=[]; % avg. area at the two conjugative nodes

%----------------User Inputs--------------------------%
e = input('Enter the total number of elements (even number): ');
n = e +1 ;      % No. of nodes = No. of elements + 1
k1 = input('Enter the value of k of Al (in W/mC): ');  
k2 = input('Enter the value of k of Cu (in W/mC): ');
q1 =input('Enter the value of heat flux entering the rod (in W/m2): ');
T_end = input('Enter the value of temperature at right end of the rod (in C): ');
d1 = input('Enter the value of diameter at left end of the rod(in mm): ')/1000;
d2 = input('Enter the value of diameter at right end of the rod(in mm): ')/1000;
L = input('Enter the value of length of the rod in m: ');

x = linspace(0,L,n);      %discretization of the domain in equal length elements

d = d1-(d1-d2)*x/L;       %diameter as a function of distance from end 1 - d is an array
A=(pi*d.^2)/4;            %area array at differental nodal points
e_l=L/e;                  %length of an element

% Building and aseembling the global conductance matrix [K]

K= zeros(n,n);
for i=1:e
    a(i) = 0.5*(A(i)+A(i+1));     % area of the element as per the described scheme
    if i<=e/2
      k=(k1*a(i)/e_l)*[1 -1;-1 1];     %element conductance matrix for Al portion of the rod
    else
      k=(k2*a(i)/e_l)*[1 -1;-1 1];      % element conductance matrix for Cu portion of the rod
    end
    % Assembly into the global conductance matrix
    K(i,i)=K(i,i)+ k(1,1);
    K(i,i+1) = K(i,i+1)+k(1,2);
    K(i+1,i)=K(i+1,i)+ k(2,1);
    K(i+1,i+1) = K(i+1,i+1)+k(2,2);
end;    

% Building and aseembling the global heat flux vector {F}
F = zeros(n,1);
Q_in = q1*pi*d1^2/4;      %total "actual" heat coming in to the system
F(1) = Q_in;                  %FEM heat coming in
q1_FEM = Q_in/a(1);

%--------Active components of set of equations--------   
%active conducatnce
K_act = K(1:n-1,1:n-1); % removing last row and col

%active force vector
F_intermediate = F(1:n-1);
F_act = F_intermediate - T_end*K(1:n-1,n);
%-------Solution------------
T=inv(K_act)*F_act;% get T profile in degree Centigrate

T_all = T(:,1);
T_all(n,1) = T_end; % Final T array

%Heat out
Q_out=K(n,1:n)*T_all;    %total "actual" heat going out of the system
q_n=Q_out/a(e);

%heat flux through each element
q = zeros(e,1);
for i = 1:e
    if i<=e/2
      q(i)=k1*(T_all(i)-T_all(i+1))/e_l;     %element conductance matrix for Al portion of the rod
    else
      q(i)=k2*(T_all(i)-T_all(i+1))/e_l;      % element conductance matrix for Cu portion of the rod
    end
end

plot(x,T_all,'ro-');
xlabel('Distance (m)');
ylabel('Temperature (Deg. C)');
title('Temperature Distribution');