% eqheat solves the heat equation using finite differences
%
% An insulated iron rod is heated to 100 C and then the ends 
%  are placed in ice baths (so the ends are fixed at 0 C).  
%  This program solves and plots the temperature distribution 
%  of the rod as a function of time.

% Based on eqheat.c from Landau/Paez and dftcs.m from Garcia.
%  Uses the "Forward Time Centered Space" (FTCS) scheme.
%  
% Programmer:  Dick Furnstahl  furnstahl.1@osu.edu
%
% Revision history:
%  26-Mar-2005 --- original version
%
clear; help eqheat;   % Clear memory and print the header

length = 100;            % length of pipe in cm
conductivity = 0.17;     % units: (cals/ s cm degree C)
specific_heat = 0.105;   % units: (cal/ g degree C)
rho = 7.87;              % density in (g/cm^3)

size = 151;                 % # of grid points in x
delta_x = length/(size-1);  % the grid spacing in x (cm)
num_steps = 20000;          % number of timesteps
delta_t = 1;                % spacing in time (s)

constant = conductivity/(specific_heat*rho) ...
                * delta_t/(delta_x^2);
if (constant < 0.5)
    disp('Solution is expected to be stable');
else
    disp('WARNING: Solution is expected to be unstable');
end

% Set initial and boundary conditions for the temperature
temp = 100.*ones(size,1);   % set temperature array to 100 C everywhere
temp(1) = 0.;      % temperature at first grid point is always 0 C
temp(size) = 0.;   % temperature at last grid point is always 0 C

% Plot variables
temp_plot(:,1) = temp(:);   % initialize the temperature plot
t_plot(1) = 0.;             % first element of time array is 0
iplot = 2;                  % initialize counter for subsequent plots
nplots = 50;                % number of snapshots to take
plot_step = num_steps/nplots; % time steps between plots
x_plot = (0:size-1)*delta_x;  % array of x-positions in delta_x steps

% Loop over the time steps
for istep = 1:num_steps
   % find new temperature from old using FTCS scheme
   temp(2:(size-1)) = temp(2:(size-1)) + ...
    constant*(temp(3:size) + temp(1:size-2) - 2*temp(2:size-1));

   % record the temperature periodically
   if ( rem(istep,plot_step) < 1 )      % every plot_step steps
       temp_plot(:,iplot) = temp(:);    % record temp(i)
       t_plot(iplot) = istep * delta_t; % set current time
       iplot = iplot + 1;               % increment the plot label
   end
end

% Plot temperature versus x and t as a wire mesh
figure(1); clf;                 % set the figure to be 1 and clear it
mesh(t_plot,x_plot,temp_plot);  % Wire-mesh surface plot
xlabel('Time');  ylabel('x');  zlabel('T(x,t)');
title('Temperature diffusion in a solid rod')
pause(1);

% Plot temperature versus x and t as a colored contour plot
figure(2); clf;              % set the figure to be 2 and clear it
contourLevels = 0:10:100;    % put contours every 10 C 
contourf(t_plot,x_plot,temp_plot,contourLevels); % contours
%cs = contourf(t_plot,x_plot,temp_plot,contourLevels); % contours
%clabel(cs,contourLevels);   % add labels to the contours
colorbar;                    % add a legend for the colors
xlabel('Time');  ylabel('x');
title('Temperature contour plot for rod');

% Plot temperature versus x and t as a smooth pseudocolor plot
figure(3); clf;              % set the figure to be 3 and clear it
cs = pcolor(t_plot,x_plot,temp_plot); % color map
shading interp;                       % remove the grid lines
colorbar;                             % add a legend for the colors
xlabel('Time');  ylabel('x');
title('Temperature plot for rod');