clc;
clear all;

%--------------Input the known------------------------%

k_Al= 200; %input('Thermal conductivity of Al (in W/mC):');
k_Cu= 389; %input('Thermal conductivity of Al (in W/mC):');  
L = input('Length of the rod (in m): ');
d=input('Diameter of the rods (in mm): ');
n=input('Number of elements in each segment of the rod: ');   %equal number of elements in each portion of the bimetallic rod               


q_in = input('Heat flux in at the left-most node (W/m^2) :');
T_node = input('Temperature of the right-most node (in Deg. C): ');

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

Le = L/(2*n);              %length of an element
Ae=(pi/4)*(d/1000)^2;      %cross sectional area of an element

K_Al=(k_Al*Ae/Le)*[1 -1;-1 1];   %elemental condutance matrix for Al
K_Cu=(k_Cu*Ae/Le)*[1 -1;-1 1];   %elemental condutance matrix for Cu

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

%Contribution from Al side
for i=1:n
    K([i:i+1],[i:i+1]) =  K([i:i+1],[i:i+1])+K_Al;
end

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

%-----------gradient matrix------------------------%
Fg=zeros(2*n+1,1);
Fg(1) = q_in*Ae;        %given boundary condition

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

T = [T_act([1:2*n]);T_node];                      %Full temperature 
q_out = single(-K(2*n+1,[1:2*n+1])*T/Ae);         %- sign as the heat is going out

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

%plotting temperatures

s=0:0.1:1;                            %dimensionless co-ordinate
for i = 1:2*n                         %evaluation for element number
    for j=1:11                        %evaluate at each dimensionless co-ordinate
        x(10*(i-1)+j)=(i-1)*Le+Le*s(j);
        Temp(10*(i-1)+j)=T(i)*(1-s(j))+T(i+1)*s(j);
    end
end

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