% Study of heat transfer and temperature of a 1x1 metal plate
%heat is dissipated through the sides and temp at infinity is t
% n(n-1) points in consideration, 
%Temperature at top end is 500*sin(((i-1)*pi)/(n-1)
%A*Temp=U
%finite difference method
% author Sidhartha Agrawal 
% BITS Pilani Goa campus India


%%
clear;
n=30;
b=525/(n*215);
t=28 ;         % temperature at infinity
% temp at top surface = 500*sin(((i-1)*pi)/(n-1))
%%
%Matrix A

%part A.1 for diagonal elements

for i=1:n*(n-1)
    if ((mod(i,n)==0) || (mod(i,n)==1) ||((i>n*n-2*n)&&(i<n*n-n)))
        if((i==n*n-2*n+1)||(i==n*(n-1)))
            A(i,i)=-1-b;
        else
            A(i,i)=-2-b;
        end
    else
            A(i,i)=-4;
    end
end



%part A.2 for other elements

for i=1:n*n-2*n
    if ((mod(i,n)==0) || (mod(i,n)==1))
        A(i,i+n)=0.5;
    else
        A(i,i+n)=1;
    end
end


% part A.3 for still other elements
for i=1:n*n-n-1
    if(mod(i,n)==0)
        A(i,i+1)=0;
    else
        A(i,i+1)=1;
    end
end

for i=n*n-2*n+1:n*n-n-1
    A(i,i+1)=0.5;
end
    
 % part A.4 for replication along the diagonal
    
for r=1:n*n-n
   for j=1:r
   A(r,j)=A(j,r);
   end
end


%%
% part C for making the constant matrix

for i=1:n*n-n
    if ((mod(i,n)==0) || (mod(i,n)==1) ||((i>n*n-2*n-1)&&(i<n*n-n+1)))
        C(i)=-b*t;
    end
end

for i=1:n
    if ((i==1) || (i==n))
        C(i)=-b*t -500*sin(((i-1)*pi)/(n-1));
    else
        C(i)=-500*sin(((i-1)*pi)/(n-1));
    end
end

%%
%solving the equations
I=inv(A);
Temp=I*C';
%%
%T=matrix containing the temperatures in terms of cordinates
for i=1:n-1
    for j=1:n
T(j,i)=Temp(j+(i-1)*(n));
    end
end

%%
% plot Surface temperature on Z-axis
%X and Y axis as cordinates of the  plate
[i,j]=meshgrid(1:1:n-1,1:1:n);
surfc(i,j,T) , xlabel('j'),ylabel('i'),zlabel('Temperature')
colormap(Hot)