% Solve Heat Eqn with Fourier Series and plot temperature
% distribution as a function of time.  -A. Ruina
% Both bar ends are at T=0
% u is the temperature.
% Choice of 3 starting temp distributions.

% The only tricky Matlab is the 'erasemode', 'xor', drawnow stuff
% described in Pratap book.

%Set up basic constants
fprintf('Calculating  ...  ') 

t_tot   =.1;    % total time of integration
ntimes=500;     % temporal resolution
numbpts = 100;  % spatial  resolution
nterms  = 100;  % number of terms in Fourier series

L=1;         %length of bar
beta=1;      %thermal diffusivity

x=(0:numbpts)*L/numbpts;  %value of x at discrete points

%Contribution to u comes from the Fourier terms
%The matrix A has one row for each term, each column is a 
%value of x
A=zeros(numbpts+1,nterms);  % initialization
for n=1:nterms
	A(:,n)=[sin(n*pi*x/L)]';
end

%Set initial temperature distribution
%fx is the initial temp distribution
%Pick from one of the three below
 %fx=x*0 +1;                 % Constant initial temperature
 fx= exp(-100*(x-.5).^2);    % Delta function, approximately
 %fx=sin(1*pi*x/L);          % sine wave 

%Find Fourier coefficients
%Crude numerical integration to find Fourier coefficients
for n = 1:nterms;
    a_n = 2 * fx*A/numbpts;	
end


u0 = A*a_n';     %Add the rows of A to get initial temp dist.
h1 = plot(x,u0 ,'erasemode','xor');  %fancy way to plot
xlabel('x,   position on bar');
ylabel('u,   temperature')
title('Temperature  vs   position, evolving in time.')

%Plot temp at a sequence of times
for t=linspace(0,t_tot,ntimes);
  %Add the rows of A, attenuated by exponential decay, to get u
  u=	A*[a_n.* exp(-beta*(1:nterms).^2 *pi^2*t/L^2)]';
  set(h1,'ydata', u);drawnow;   % fancy plot, continued	
end

fprintf(' done.\n\n')