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/3/Le)*[7 -8 1;-8 16 -8;1 -8 7];   %elemental condutance matrix for Al
K_Cu=(k_Cu*Ae/3/Le)*[7 -8 1;-8 16 -8;1 -8 7];   %elemental condutance matrix for Cu

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

%Contribution from Al 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_Al;
end

%Contribution from Cu side
for i=n+1:2*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_Cu;
end

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

%-------------solution------------------------%
K_act = K([1:4*n],[1:4*n]);                       %as the last node temperature is given
Fg_act = Fg([1:4*n])-K([1:4*n],4*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:4*n]);T_node];                      %Full temperature 
q_out = single(-K(4*n+1,[1:4*n+1])*T/Ae);         %- sign as the heat is going out

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

%plotting temperatures

s=0:0.05:1;                            %dimensionless co-ordinate
for i = 1:2*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);
xlabel('x (m)');
ylabel('Temperature (Deg C)');