%  ============================================================
%  chimney.m -- Compute and displace profiles of temperature in 
%               chimney cross section.
%  ============================================================

%  Setup working array and boundary conditions along
%  internal/external walls.

 T = zeros(9,9);

 for i = 5:9;
     T(i,5) = 200;
 end;
 for i = 1:5;
     T(5,i) = 200;
 end;
 for i = 6:9;
 for j = 1:4;
     T(i,j) = NaN;
 end;
 end;

% Loop over internal nodes and compute new temperatures
           
 counter   = 0;
 maxchange = 200;
    while (maxchange > 1)
        counter = counter+1;
        maxchange = 0;
        k=5;
        l=4;
        for c = 6:8;
            newtemp    = 0.25*(2*T(l+1,c)+2*T(l,c+1));
            tempchange = newtemp - T(l,c);
            maxchange  = max(maxchange,abs(tempchange));
            T(l,c)=newtemp;
           for r=k:8
               newtemp    = 0.25*(T(r,c-1)+T(r,c+1)+T(r-1,c)+T(r+1,c));
               tempchange = newtemp - T(r,c);
               maxchange  = max(maxchange,abs(tempchange));
               T(r,c)= newtemp;
           end
           newtemp    = 0.25*(T(9,c-1)+T(9,c+1)+2*T(8,c));
           tempchange = newtemp - T(9,c);
           maxchange  = max(maxchange,abs(tempchange));
           T(9,c)=newtemp;
           l=l-1;
           k=k-1;
        end
        counter;
        maxchange;       % to view counter or maxchange remove 
   end
     
% Compute reflected temperature 

  for i = 2:4
      for j = 1: 11-i
          T(i,j) = T(10-j,10-i);
      end
  end

% Print temperature array.

T

% Plot temperature contours.

contour(T) 
hold;

% Now overlay perimeter of chimney section on contours.

perim = [ 1 , 1;
          9 , 1;
          9 , 9;
          5 , 9;
          5 , 5;
          1 , 5;
          1 , 1 ];
     
plot(perim(:,1),perim(:,2),'w');
text(1.1,5.3,'Temp = 200.0');
text(7.0,1.3,'Temp = 0.0');

%  ============================================================
%  The End!